Relationships between substrate , surface characteristics , and vegetation in an initial ecosystem

We investigated surface and vegetation dynamics in the artificial initial ecosystem “Chicken Creek” (Lusatia, Germany) in the years 2006–2011 across a wide spectrum of empirical data. We scrutinized three overarching hypotheses concerning (1) the relations between initial geomorphological and substrate characteristics with surface structure and terrain properties, (2) the effects of the latter on the occurrence of grouped plant species, and (3) vegetation density effects on terrain surface change. Our data comprise and conflate annual vegetation monitoring results, biennial terrestrial laser scans (starting in 2008), annual groundwater levels, and initially measured soil characteristics. The empirical evidence mostly confirms the hypotheses, revealing statistically significant relations for several goal variables: (1) the surface structure properties, local rill density, local relief energy and terrain surface height change; (2) the cover of different plant groups (annual, herbaceous, grass-like, woody, Fabaceae), and local vegetation height; and (3) terrain surface height change showed significant time-dependent relations with a variable that proxies local plant biomass. Additionally, period specific effects (like a calendar-year optimum effect for the occurrence of Fabaceae) were proven. Further and beyond the hypotheses, our findings on the spatiotemporal dynamics during the system’s early development grasp processes which generally mark the transition from a geo-hydro-system towards a bio-geo-hydro system (weakening geomorphology effects on substrate surface dynamics, while vegetation effects intensify with time), where pure geomorphology or substrate feedbacks are changing into vegetation–substrate feedback processes.


Introduction
While a lot of studies on ecosystem development have been conducted in mature ecosystems (e.g.Campbell et al., 2007;Ellenberg et al., 1986;Fränzle et al., 2008;Pennisi, 2010) our information about initial systems is comparably weak.This is remarkable because initial ecosystems are usually less complex compared to mature systems.Therefore, the study of the development of initial ecosystems could be very helpful to achieve a better understanding of the complex relationships and feedback mechanisms which typically arise over time (Jørgensen et al., 2000).Tracing the development of young ecosystems and observing how new relationships and feedbacks emerge with increasing complexity would help to get a better insight on key processes and a basic understanding of their interactions.
Initial ecosystems are also important from a very practical point of view.Worldwide, human activities regularly leave behind huge bare areas (Walker and Willig, 1999;Bradshaw, 1983;Hüttl and Weber, 2001;Schaaf, 2001;Zikeli et al., 2002), which create serious environmental challenges due to Published by Copernicus Publications on behalf of the European Geosciences Union.P. Biber et al.: Relationships between substrate, surface characteristics, and vegetation their ecological instability and low productivity.Here, further knowledge is urgently needed for finding optimal ways to transform such landscapes into landscapes that can be used in a sustainable way by society.
The land surface is an interface where geomorphological and biological processes connect.It can also be seen as an interface where pedosphere, biosphere, hydrosphere and atmosphere are strongly interlinked (Brantley et al., 2007).During early developmental stages, the evolution of a rill system which channels surface runoff and material transport probably is the most visible outcome of such interactions.At the very beginning of ecosystem development, starting with bare ground, these interactions are expected to be especially important (Schaaf et al., 2012).
Scientifically investigated areas, where the first interactions between atmosphere, pedosphere and biosphere are just starting to develop, are rare.Among them are areas created by volcanic activity (Bishop, 2002;Dahlgren et al., 1999;del Moral and Wood, 1993;Müller-Dombois and Fosberg, 1998;Friðriksson, 2005) or areas that become exposed after glacier retreat (Cooper, 1923;Matthews, 1992).However, information about the initial conditions at "point zero" is generally incomplete.
The creation and examination of artificial areas is one approach to close the existing knowledge gap on processes determining early ecosystem development.Ideally, such areas are complete water catchments with a well-defined and if possible homogeneous substrate, and with negligible traces of life and therefore without any successional history.Their value for understanding the emerging interactions between atmosphere, pedosphere, hydrosphere and biosphere cannot be overestimated (Schaaf et al., 2011).The 6 ha artificial catchment "Chicken Creek", established in 2004/05 in the open-cast mining area Welzow-Süd near Cottbus, eastern Germany, was designed to represent the theoretical ideal conditions as close as possible.Here, over a period of seven years, so far, a broad range of key properties of geomorphology, vegetation succession and ecosystem development have been studied with high spatio-temporal resolution.
Combining interdisciplinary data from the Chicken-Creek catchment, our goal within this study is to capture the developing complexity in an initial ecosystem by scrutinizing the following overarching hypotheses.
-Initial substrate characteristics determine the terrain surface's structure (rills, local relief energy) and its change.The degree of determination of these processes weakens with the ecosystem's development (H1).
-Surface and substrate characteristics determine the plant species groups to be found initially (H2).
-Increasing density of vegetation reduces the erosion and therefore the degree of change of terrain surface structure (H3).
More specifically, in the context of H1 we expect to find the occurrence of rills, the local relief energy, and the terrain surface height change from one year to another being timedependently correlated to the annual groundwater level, the near-surface substrate grain size, and other substrate characteristics as measured after catchment construction.
From the perspective of H2, vascular plant species are grouped by ecological and morphological criteria into annual, herbaceous, grass-like, and woody, and affiliation to the N-fixing Fabaceae family.The occurrence of a given plant species group is expected to depend on the same variables as mentioned above.Moreover, H2 comprises total plant cover, overall vegetation height and density as goal variables which are assumed to be related to the above-mentioned explanatory variables in terms of a progressive development as well.
In contrast to H1, where only groundwater and substrate characteristics are related to the change of surface height, H3 includes a measure of vegetation density as an explanatory variable.

The artificial catchment Chicken Creek
The artificial catchment, Chicken Creek, with an area of 6 ha was constructed in the open-cast mining area of Lusatia, Germany (51.6049 • N, 14.2667 • E) in [2004][2005].It is a 2-4 m layer of post-glacial sandy to loamy sediments overlying a 1-2 m layer of Tertiary clay which forms a shallow pan and seals the whole catchment at the base.No further measures of restoration like planting, amelioration or fertilization were carried out; natural succession and undisturbed development is allowed.Due to the artificial construction, subsoil boundary conditions of this site are clearly defined including well documented inner structures as compared to natural catchments.
The region around the catchment is characterized by a temperate climate with a sub-continental character and comparatively low precipitation.During the last five decades, the mean annual temperature amounted to, on average, 8.9 • C (January mean: −0.8 • C, July mean: 18.4 • C), the annual precipitation to 563 mm (Gerwin et al., 2009).
Up to now, the set of plant species that has been observed during succession is very similar to what field botanists can find in the region (e.g.Felinks, 2000).Propagules of course are brought in by seed rain, but there has also been a small but undeniable seed bank (Zaplata et al., 2010(Zaplata et al., , 2011a)).
In the initial situation, organic substance was not completely absent from the soil, concentrations however, were very low (about 2 mg g −1 substrate).Radio-carbon dating yielded ages between 3000 and 16 000 which indicate the organic substance being mostly of fossil origin (Gerwin et al., 2009;Risse-Buhl et al., 2013).
For more details on the construction process, initial site conditions as well as the monitoring program carried out since 2005 see Gerwin et al. (2009Gerwin et al. ( , 2010Gerwin et al. ( , 2011) ) and Schaaf et al. (2011Schaaf et al. ( , 2012)).

Data used in this study
A 20 m × 20 m grid for sampling and orientation purposes was established across the catchment area (Fig. 1).The data used in this study cover the sub-area between the grid-point lines A and Q as shown in Fig. 1 which amounts to approximately 5 ha.The other part, the area around the pond in the lower part of the catchment, was excluded from our study since it is influenced by a high, constant water table and evolved into a semi-aquatic ecosystem with different site conditions.
For the statistical scrutiny of our hypotheses a combined data set comprising terrain surface information from terrestrial laser scanning (TLS), vegetation information from field records, soil and groundwater information was used.The data sources have different spatial and temporal resolutions.Soil samples were taken only once at the beginning of the study, vegetation records were made once every year in summertime and TLS was done at seven different times between 2008 and 2011.
The above-mentioned 20 m × 20 m grid points form the basic units of each analysis carried out.A spatial data fusion was thus done in order to assign the different data sources to the grid (see below).Temporally, vegetation record data were assigned to the same year's laser scan data.Thus, dependent on the analysis of interest, the same vegetation record data were assigned to both laser scan measurements of a year.Groundwater measurements were assigned as annual means to each grid point.As there was only one survey of soil properties, a specific temporal assignment was not possible.Thus, for interpretation purposes they represent initial conditions and not a temporal development.
As the catchment is sloped downwards from row A to row Q, the distance of any grid point to row A is indicated as DISTA.The time since January 2005 is encoded in the variable TIME with month as yearly fractions (e.g.6.42 for May 2011) and CALYEAR as the calendar year.In the following text, the resulting time component of the variables used in this study is identified by the index j and the spatial location by the index i.In Table 1 we list all variables and indexes together with their definitions, while in Table 2 we present their summary statistics (more details below).All response variables were chosen and calculated a priori, based on considerations about ecological interpretability and plausibility.

Vegetation
Altogether 119 permanent plots for long-term vegetation monitoring were established at the above-mentioned 20 m × 20 m grid points (Zaplata et al., 2010;Fig. 1).Each plot had a square shape of 5 m × 5 m with the respective grid www.biogeosciences.net/10/8283/2013/Biogeosciences, 10, 8283-8303, 2013 point in its centre.At each plot all vascular plant species were recorded including cover estimates for each species in 30 distinct percentage classes (Zaplata et al., 2013).Bryophyte covers were estimated in the same way, distinguishing mosses and liverworts.For this study, a time series with 6 yr of annual monitoring (2006)(2007)(2008)(2009)(2010)(2011) is available.The plant species were grouped by (i) lifespan according to Rothmaler (2000), but with only the two categories "annual" and "perennial"; (ii) life forms "grass-like", "herb" or "wood" according to Rothmaler (2000), species of the genus Rubus which in general form woody stems were thus labelled as "wood"; and (iii) affiliation to the Fabaceae family ("Fabaceae" versus "no Fabaceae").Fabaceae species are of special interest, because they are the only N-fixing vascular plants found on the catchment.Other N-fixing organisms are cyanobacteria, which are part of biological soil crusts, also found at the study site (Fischer et al., 2012).However, in this context Fabaceae is the most important group by far.See Table A.1 in the Appendix for each vascular plant species' group affiliation.
In Table 3 we summarize the covers separately for the above-mentioned groups and for all plants together.These group-wise and overall covers were obtained by summing up the single species' covers as recorded in the field.The total plant cover increases from an average of 2.0 % in 2006 up to 105.7 % in 2011.The latter value, being greater than 100 %, indicates an overlap of the areas covered by different species.This dramatic rise occurs rather homogeneously across the area, albeit with a slightly higher cover on its central longitudinal axis (Supplement Fig. S1).While from 2006 to 2008 the vascular plants (annuals and perennials) almost exclusively contribute to the total cover, the further increase of total cover is carried by a massive increase of the bryophytes between 2010 and 2011 (Elmer et al., 2013).The vascular plants' cover, however, seems to stabilize between 30 and 40 % during 2009-2011.
The annual plants reach a distinct maximum cover of more than 30 % in 2009 which considerably decreases in the following years.A similar development, but with a less steep decline, can be seen for the herbaceous plants because this group largely overlaps with the annuals (Supplement Figs.S2, S3).Since 2008, Fabaceae occur with very similar cover values as the annuals do.This is due to the annual species Trifolium arvense L. (Fabaceae) which is the most prominent annual species from 2008 to 2011 (cf.Zaplata et al., 2011b, Supplement Fig. S4).
The perennial plants exhibit a steady area-wide increase of cover during the whole observation period, with however still considerably lower covers than for the annuals in 2011 (Supplement Fig. S5).On a lower level, the grass-like plants show a mostly parallel development, seeming to stabilize in 2010 and 2011, with Supplement Fig. S6 suggesting strong increases at the southeast edge (near the pond) and a slight decrease across the rest of the area.Compared to the other plant groups, the covers of woody plants seem negligible so far.However, their average cover increased exponentially and their presence on the whole area (Supplement Fig. S7) might be a precursor of their future dominance.
In order to capture the different groups' relative dominance and to eliminate the time trend which is inherent in the absolute cover values, we related their covers for each grid point and each observation to the respective total plant cover COVTOT for further analysis.This resulted in the proportions of cover for annual plants (PROPANN), grasslike plants (PROPGRASS), woody plants (PROPWOOD) and Fabaceae family (PROPFAB) (cf.Tables 1, 2).As Table 2 shows, the herbaceous plants are the dominating vascular plants throughout the observation period, however with a pronounced decrease of their share from between 0.7 and 0.8 down to about 0.4 and 0.3 in 2010 and 2011.Annual plants follow the same pattern, albeit on a slightly lower level (cf.Supplement Figs.S8, S9).Grass-like plants increase their share from 0.2 to almost 0.3 in 2008.After that, they also drop down, but by far not as dramatically as observed for the annuals and herbaceous plants (cf.Supplement Fig. S10).The Fabaceae, almost not present in 2006, extend their share, starting from the area's western part (Supplement Fig. S11), up to nearly 0.6 in 2009 and drop -in line with the annuals -back to 0.2 in 2011.The low shares of all these groups in 2010 and 2011 come from the above-mentioned vast increase of bryophyte covers.Only the woody plants' share of the total cover keeps steadily increasing, albeit on a very low level so far (cf.Supplement Fig. S12).

Soil and groundwater
After completion of the catchment construction, soil samples were taken at each grid point at a depth of 0-30 cm between October 2005 and April 2006 (see Fig. 1).The samples were air-dried, passed through a 2 mm mesh and analysed for pH (water extract, Beckmann pH34 glass electrode and WTW pH537), electrical conductivity (EC, Hanna HI 8733 and WTW LF537), texture (sieving and sedimentation procedure with Köhn pipette method), total content of carbon (C T ), nitrogen (N T ) and sulfur (S T , elementar analyser Vario EL III), inorganic (carbonate, Scheibler calcimeter) and orwww.biogeosciences.net/10/8283/2013/Biogeosciences, 10, 8283-8303, 2013 ganic carbon (C org , calculated as difference C T -CaCO 3 -C) content (Sparks, 1996;Dane and Topp, 2002).We assumed the soil data from a grid point being valid for the whole vegetation square.Variables that turned out to be promising for this study from a theoretical point of view as well as in visual exploratory data analyses were the percentage of C org (CORGP), the percentage of gravel content (GRAVCONT), and MSAND, the percentage of medium sand (see Tables 1, 2).Hereby, the gravel fraction of a soil is defined by grain sizes of more than 2 mm.C org -mostly from fossil sources -was present in the initial situation, albeit with low concentrations only (Risse-Buhl et al., 2013).These generally low concentrations of CORGP show considerable variation between 0 and almost 0.6 % (Table 2).Spatially, CORGP is not evenly distributed, higher concentrations are found in the south-west half of the area (Supplement Fig. S13).Mean gravel content (GRAVCONT) is about 15 % (Table 2) with a slight tendency towards higher values again in the south-west zone (Supplement Fig. S13).At an average of 47 % MSAND shows a homogeneous spatial distribution across the catchment (Table 2, Supplement Fig. S13).
A total of 21 groundwater gauges were initially installed across the catchment.Nine of them were equipped with pressure transducers to register groundwater levels automatically.The levels at the other gauges were measured manually every two weeks.During the following years, five additional groundwater gauges were installed.The positions of all 26 gauges can be taken from Fig. 1.See Biemelt et al. ( 2010) and Supplement Fig. S14 for more details.We simply attributed the mean annual distance of the groundwater level from the surface (variable GAUGE, see Table 1) from the nearest measurement shaft to each grid point, with an average distance of 19 m and maximum distance of 53 m from a square midpoint.We used that simple method because no complex assumptions are made as is the case for any interpolation method.
Since 2006, GAUGE shows a declining trend from on average 1.8 down to 0.9 m in 2011, indicating rising groundwater levels over time.Their spatial distribution (Supplement Fig. S14) is homogeneous in the beginning with higher levels only in the lowest part (south-east) with a slight tendency towards drier conditions at the south-west edge developing during the observation period.

Terrestrial laser scans
We used a terrestrial laser scanner (mod.Riegl LMS Z420i) in last target mode to measure 3-D ground surface and vegetation height and density simultaneously.In order to keep the impact on the catchment at a minimum, we only measured from 13 permanently fixed scan positions (Fig. 1), each geo-referenced with 30 DGPS-fixed points (DGPS, differential global positioning system), which allowed us to keep the standard deviation of the geo-referencing error below 2.5 cm.
The scan positions were spatially arranged in a way that allowed us to maintain a horizontal measurement resolution of at least 10 cm × 10 cm at a hypothetical horizontal ground surface.Given the limited number of scan positions, this was only possible by mounting the laser scanner on a portable 6 m tower.This height is sufficient for achieving the desired minimum resolution but not too high in order to hit existing vegetation mainly from more lateral and not so much from vertical directions.This is important, because vertically oriented vegetation like grasses is more likely detected from the side than from above and the side view is better suited to detect vegetation layers under emerging woody plants.See the method comparison by Schneider et al. (2012) for more technical details.
We conducted seven area-wide laser scan campaigns between 2008and 2011(September 2008, Mai 2009, August 2009, April 2010, October 2010, February 2011, and May 2011).Although it was not easy to find appropriate time windows (for such a campaign several days with clear weather and low wind are required) the measurements were timed so that at least in the years 2009 and 2010 roughly the catchment's state at the beginning and the end of the growing season could be covered.With our scans we covered the whole catchment amounting to 6 ha altogether.
The raw laser scan data are three-dimensional point clouds, each point indicating the position where the laser beam hit an object.In order to separate vegetation from ground surface we divided the covered area into about following form: with a, b, c being parameters to be estimated and ε representing the error term.Equation (1) thus represents a plane that after fitting can be interpreted as a smoothed description of the terrain surface in the cell of interest's proximity defined by the 6 m radius as explained above.If the actual ground level z of the central cell is more than 5 cm lower than the one estimated from Eq. ( 1), then this cell is assumed to belong to a rill.Visual comparisons with aerial photographs from the catchment show a good agreement with this rill detection method.
In addition to this 6 m neighbour analysis, we used the same method with a 0.8 m radius around the cell of interest and calculated the difference between the greatest and the smallest vertical deviation from the regression plane.This can be interpreted as a description of relief energy in the nearest neighbourhood.This potential energy is used as an indicator if the cell is located on a jagged or smooth surface and thus how large the influence of surface runoff events might be or how likely a successful seed deposition may happen.
Changes of the ground surface level inform about erosion and deposition processes.For each cell we calculated the difference between the surface level at a given and a previous survey, requesting a survey-to-survey time distance closest to one year.Negative numbers, i.e. decreasing surface levels, indicate erosion while positive numbers indicate deposition.For our analysis the changes were always attributed to the first of the two surveys, in other words to the status before the change.
To each 20 m × 20 m grid point we attributed every laser scan grid cell with its midpoint inside the corresponding 5 m × 5 m vegetation square and used only those cells thereafter.This resulted in about 100 laser scan cells per vegetation square and implied two options for data evaluation.First is the extraction and evaluation of sum and mean values from the laser scan data, one number per vegetation square, such as surface roughness or percentage of rill cells.Second is an analysis on the level of single laser scan grid cells.However, preliminary analyses showed that the latter tended to be less revealing than the former while consuming disproportional computing time.Thus, we focussed on the former, more aggregated option.
For further analysis, the following variables were calculated from the laser scan data, resulting in one number per grid point and survey (cf.Table 1): VEGHEIGHT and VEG-DENS are the mean vegetation height and the mean vegetation density of all laser scan grid cells on a vegetation square.VEGHD is their product.RILLR is the proportion of rill cells as defined above on a given vegetation square while RELEN is the average relief energy of all laser scan cells on one vegetation square.DHEIGHT represents the average surface height change of all laser scan cells belonging to a given vegetation square.
On average VEGHEIGHT more than doubles from 2008 to 2011 (0.35 to 0.80 m, Table 2).The highest vegetation concentrates in the southern and south-west part of the catchment in 2011 (Supplement Fig. S15).The vegetation density, VEGDENS, on average remains constant during the observation period, albeit with changing spatial distribution (Table 2, Supplement Fig. S16).Its overall constancy despite increasing heights confirms its adequacy for expressing vegetation density independent from plant size.Accordingly VEGDH, the product of both, doubles from 2008 to 2011 in a similar way as VEGHEIGHT, being greatest in the southern corner and the east side of the area, where both, height and density are high (Table 2, Supplement Fig. S17).
The average proportion of rill cells per vegetation square, RILLR, does only slightly increase from 0.10 to 0.16 from 2008 to 2011 (Table 2).Being concentrated on the northern, western and southern edges of the investigated area from 2008 to 2010 a redistribution towards the south-eastern part of the catchment seems to have occurred between 2010 and 2011 (Supplement Fig. S18).Similar to RILLR, the local relief energy, RELEN, stays nearly constant between 0.12 and 0.13 from 2008 to 2010 but almost doubles from 2010 to 2011 accompanied with a spatial redistribution of the highest values towards the eastern edge (Table 2, Supplement Fig. S19).The surface height change, DHEIGHT, undergoes a reversion during the observation time span.In average, one vegetation square lost 3 cm of terrain height between 2008 and 2009, from 2009 to 2010 there was on average no change, and from 2010 to 2011 we observe an average surface height increase of 4 cm (Table 2).A spatial view, as provided in Supplement Fig. S20, shows that erosion first dominated throughout the catchment, being strongest at the western edge.Later we observe a central ridge where terrain height slightly increases while erosion weakens but continues along the western and the eastern edges.From 2010 to 2011 deposition dominates the whole area, being strongest in the north-western part.

Statistical evaluation
We chose generalized linear models (GLMs) and generalized additive models (GAMs) as the basic method for our statistical evaluations, because they allow us to deal with different distribution families for the response variable.In addition, GAMs offer a convenient way to incorporate nonlinear effects into the analyses.However, as we take a vegetation square as our basic observation unit, several subsequent observations are available per vegetation square.We covered the resulting autocorrelation by introducing a random effect on vegetation square level.This enlarges GLMs to generalized linear mixed models (GLMMs) and GAMs to generalized additive mixed models (GAMM's).See Zuur et al. (2009) for a detailed description of these model types and Pinheiro and Bates (2000) for further information about mixed models.Preliminary tests showed that including spa- tial correlation between vegetation squares did not improve our models.
The core of a GLMM in the context of this study is a linear predictor function η of the following form: , where X 1 . . .X q , is a set of q explanatory variables.The indices i and j denote the j th observation on the ith vegetation square.α and β are regression parameters and a i is a vegetation-square specific random effect.In case of a GAMM, the predictor function enlarges to where, in addition, K 1 . . .K n is a set of n explanatory variables, and f 1 . . .f n is a set of nonparametric smoother functions (cf.Zuur et al., 2009).
The way, that η is transformed into an expected value µ ij for the actual response variable Y ij is defined by the link function g.
If the response variable Y ij is assumed to be normally distributed (Y ij ∼ N (µ ij , σ 2 )) then the usual link function is identity: This is the first of three options we used in this study.
In cases where Y ij was proportional data with a tendency towards over-dispersion, a quasi-binomial distribution for and the dispersion parameter τ not being fixed to 1 and a logit link was used.When Y ij > 0, the assumption of a gamma distribution for Y ij with var(Y ij ) = µ 2 ij ν and ν being the dispersion parameter turned out to be useful in combination with a logarithmic link: Promising explanatory variables for each single regression analysis were pre-selected by theoretical considerations and, as proposed by Zuur et al. (2009), based on descriptive exploratory data analyses such as pair-plots and varianceinflation factors.If plausibility, visual data inspections or residual plots suggested (partly) nonlinear relationships, a GAMM was chosen instead of a GLMM.In each analysis, we started with the full set of pre-selected explanatory variables, leaving them out stepwise, by AIC comparison and refitting each time (Zuur et al., 2009).Model assumptions were confirmed by graphical displays.
Presenting the most important results of a model fit is straightforward for a GLMM (Eq.2).Presenting the linear parameters' β 1 , . .., β q estimates together with their standard errors and significance levels is standard.In the case of a GAMM, this can be done in exactly the same way for the linear components of the model.However, the components which are expressed by nonparametric smoother functions (f 1 (K 1ij ), . .., f n (K nij ) in Eq. 3) can only be displayed visually.Thus, for a GAMM, a table is needed for presenting the linear components and a diagram is needed for each component that is represented by a smoother function.
In order to scrutinize H1, "do the initial geomorphological and substrate characteristics determine the surface's structure and properties of the terrain?", we investigated several response variables (see Table 1): -RILLR: the proportion of rill cells per vegetation square and measurement.As RILLR represents proportional data, the quasibinomial distribution with logit link was used.
-RELEN: the mean value of the local relief energy of all laser scan cells belonging to one vegetation square.Assuming a gamma distribution combined with the logarithmic link matched the data well from a theoretical point of view and yielded good model fits for this response variable.
-DHEIGHT: the ground's height change (see Table 1 for details) was investigated with a regression model that assumes normal distribution of the response variable and uses the identity link.
Due to the definition of H1 we did not include any vegetation properties as explanatory variables in these analyses.
In the context of H2, asking if surface and substrate characteristics determine the plant species groups to be found initially, the following set of vegetation properties was chosen as response variables (see Table 1): -PROPANN: the portion of annual plants related to the total coverage per vegetation square and survey.
-PROPHERB: same for the portion of herbaceous plants.
-PROPFAB: same for plants belonging to the nitrogenfixing Fabaceae family.
-PROPWOOD: same for woody plant species.For all these proportions, assuming quasi-binomial distribution and using a logit link turned out to be the best option.In addition, we also used two laser-scanned vegetation properties as response variables:

Biogeosciences
-VEGHEIGHT and VEGDENS: the mean vegetation height and vegetation density per vegetation square and survey.
The former was investigated assuming normal distribution and using the identity link while for the latter quasi-binomial distribution with logit link was applied.
For testing H3, if increasing density of vegetation reduces the erosion, we used DHEIGHT: the laser-scanned surface height changes as a response variable.
We did so in the same way as for H1, but tested several ways of including the laser-scanned vegetation properties VEGHEIGHT and VEGDENS as additional explanatory variables.The other explanatory variables were those that showed significance in the context of H1.
The overall important components of temporality and spatiality are included in the variables TIME and DISTA.The final models with the final set of explanatory variables are shown in the results section.Table 1 explains all variables used in this study, and Table 2 presents their means, minima and maxima for each observation year.
For all statistical evaluations we used the free software R version 2.15.1 (R Core Team, 2012), namely the package mgcv (Wood, 2006).3 Results

Influence of geomorphological and substrate characteristics on surface structure (H1)
The final model for describing the dependency of the portion of rill cells, RILLR, was a GAMM with logit link (cf.Sect.2.3, Eqs. 3, 6): The indices i and j , and the random effect a have the same meaning as in Eq. ( 2), and E(RILLR) represents the expected value of RILLR.CORGP is the percentage of organic carbon in the soil, GAUGE is the mean annual groundwater level in m below surface, and f 1 (DISTA) is a smoother function that represents the effect of DISTA (see Table 1).All explanatory variables show highly significant (p < 0.001) effects.
Figure 2 shows the effect of DISTA.Clearly there is a general trend of more likely finding rills the more one moves further downslope, however with a minimum between 50 and 150 m from row A. Table 4 lists the parameter estimates for α and β 1 , . .., β 4 .The parameters for the variables CORGP, GAUGE, and TIME are greater than 0, indicating that the probability of encountering rills increases with the percentage of organic carbon (as measured in 2005), a lower annual groundwater level (higher GAUGE) and with time.However, there is also a significant interaction that indicates the influence of CORGP weakens with time.
Testing the local relief energy, RELEN, the right hand side of the final GAMM model was the same as in Eq. ( 8 The results as presented in Table 5 and Fig. 3 show a strong affinity between the local relief energy and the rill probability.Again, relief energy generally increases with increasing DISTA, however with a minimum at about 100 m from row A (Fig. 3).Relief energy also increases with CORGP (Table 5), lower groundwater levels (higher GAUGE) and, expectedly, increases with TIME (parameters β 1 , . . ., β 3 > 0), however the negative value obtained for β 4 shows a weakening influence of CORGP with time.
The change of the surface height, DHEIGHT, in relation to non-vegetation variables only, yielded a final model of the following GLMM structure:   6. Significances are partly weaker than in the models shown above.The effect of MSAND (β 1 < 0) means that more material (volume) erodes at places where more medium sand was found in 2005.Lower groundwater levels are connected with greater material export (β 2 < 0), the more rills we find, the more material is exported (β 3 < 0).There is a weak, but highly significant tendency towards a surface height increase -indicating deposition effects -with increasing DISTA (β 4 > 0) while time has no significant isolated effect, how-ever the analysis reveals significant time effects in interaction with other variables.Parameter β 6 , being significantly greater than 0 indicates that the above-mentioned effect of MSAND weakens with time, and the same is true for the effect of RILLR (β 7 > 0).A reverse effect is observed for DISTA in interaction with time (β 8 < 0), which describes well the spatial change of DHEIGHT as visualized in Supplement Fig. S20.

Influence of surface and substrate characteristics on vegetation structure (H2)
With the proportion of annual plants, PROPANN, as the response variable the following final model resulted (GAMM with logit-link): With β 1 < 0, Table 7 shows a negative correlation between CORGP and the proportion of annual plants.The non-linear effect of COVTOT, the respective vegetation plot's total degree of coverage in percent (Fig. 4) reveals an optimum at roughly 50 % while the increase at very high degrees of coverage is not significant (broad confidence interval).The proportion of annual plants decreases with the calendar year (Fig. 4), while the probability of finding annual plants on a vegetation plot increases with lower groundwater levels (higher GAUGE) roughly in the shape of a saturation curve (Fig. 4).

Variable
The effects observed for PROPHERB are virtually the same as obtained for annual plants (Fig. 5, Table 8), an optimum pattern for COVTOT, a decreasing curve for the calendar year and a saturation curve for the groundwater level, and a negative correlation with CORGP.
The scrutiny of the proportion of grasses (PROGRASS) yielded an even simpler model: Compared to the herbaceous plants, the grasses show an opposing trend.Even though the effect of CORG is just not significant (p = 0.0633, Table 9), the grasses seem to be connected to places with higher organic carbon content in the soil.As Fig. 6 shows, the coverage of grasses decreases strongly with the total degree of coverage (COVTOT), while it increases with the calendar year.This indicates that grasslike plants preferentially establish where the total cover is not too high.
For the proportion of woody plants (PROPWOOD) a simple GLMM resulted: Only two explanatory variables show a significant relationship in this context.Trivially, the proportion of woody 57 1102 Figure 6 1103 Fig. 6.Nonlinear effects of the total degree of coverage (COV-TOT), and calendar year (CALYEAR), on the proportion of grasses (PROPGRASS) (see Eq. 13).Shaded: 95 %-confidence area.
plants increases with time (CALYEAR).Less expectedly, the probability of finding woody plants decreases as the total coverage on a vegetation square increases (Table 10).
The final model for the proportion of Fabaceae (PROP-FAB) resulted in the following structure: with GRAVCONT representing the soil's gravel content in percent as measured in 2005.
As Table 11 shows, there is a positive correlation between lower groundwater levels and the occurrence of Fabaceae (β 1 > 0).Up to gravel contents of 10 % and up to total degrees of coverage of about 20 %, their proportion is positively correlated with both variables (Fig. 7).Remarkably, the occurrence of Fabaceae over time showed a maximum in the years 2008 and 2009 (Fig. 7).VEGHEIGHT, the mean vegetation height, as obtained from laser scans was described with the following model: Figure 8 (left) shows that vegetation height increases with time, expectedly in an exponential pattern as it is typical for initial growth processes.The same figure (Fig. 8, right) reveals a clear trend of increasing vegetation height with the distance from grid-row A. Depth to groundwater table and coverage degrees are positively correlated with vegetation height (β 1 > 0, β 2 > 0, Table 12).No significant relations of other variables with the laserscanned vegetation density, VEGDENS, could be identified.

Influence of vegetation structure on surface structure (H3)
In order to scrutinize H3 we took the model formulation from Eq. ( 8) which relates the change of the laser-scanned DHEIGHT, with soil and surface variables, but added different combinations of the laser-scanned VEGHEIGHT, VEG-DENS, and TIME, all referring to the beginning of the period during which DHEIGHT occurred.We found that including the product of VEGHEIGHT * VEGDENS = VEGDH in interaction with time yielded the best model: Compared to the model without vegetation data as explanatory variables (Eq. 10, Table 6), the parameter estimates for β 1 , . .., β 8 are very similar (Table 13), thus, there is no contradiction between this and the other model.However, the nonlinear vegetation effect pronouncedly shows a vegetation impact that counteracts erosion.The higher and denser the vegetation (i.e. the greater VEGDH), the stronger this effect, which also changes with time (Fig. 9).In the first observation period (starting September 2008) the correlation is somewhat weak and almost linear but in the following two periods, the curve becomes more and more bent upwards, the VEGDH value where the bend starts moving down to smaller VEGDH values from the second to the third observation period.This indicates that even with the same VEGDH values the stabilizing effect of vegetation increases with time.

Discussion
Surface structure in the study area is mainly influenced by geomorphological characteristics.The slope degree of the area is not the same throughout the area, but shows a maximum in the lower part (rows L to P in Fig. 1).This is reflected (or at least indicated) by the pronounced and partly nonlinear effects of the downhill distance from grid-row A (DISTA) on the frequency of rills (RILLR), the local relief energy (RELEN), and the ground surface height change (DHEIGHT).The pronounced minima for the downhill distance effect (DISTA) for the rill frequency and the local relief energy (Figs. 2, 3) seem to indicate a special local situation: they are located just downhill an area where, due to the construction works of the catchment, many rills take a perpendicular course to the main slope direction (cf.Schaaf et al., 2013, Supplement Figs. S21, 22, 23).
Gullies channelling surface runoff results in a further incision both in depth and uphill -a positive feedback, which is strongly supported by our model for the ground surface height change (DHEIGHT, Table 6), especially by the signif-icant effects of the rill frequency (RILLR) and the interaction effect of the downhill position (DISTA) with time.
Although the portion of organic carbon (CORGP) is strongly correlated with rill intensity, it seems most likely that this parameter is a proxy for the general substrate properties.The same could be true for parameter MSAND.According to Gerwin et al. (2010), CORGP strongly correlates with clay content.In general, soil texture in the western part of the catchment is dominated by loamy sands, whereas in the eastern part pure sands dominate.These differences are reflected in the spatial extent and form of the erosion gullies being narrower and deeper in the western part, but shallower and wider in the eastern part (Schaaf et al., 2012).Thus, in general, our results indicate positive feedback processes governing rill formation whereby, however, different rill types evolve on substrates with different properties.
Surface runoff dominated catchment hydrology at least in the early years, resulting in gully formation and transport of large sediment amounts that were deposited in the lower part of the hill slope and in the pond (Schaaf et al., 2012).While doing so surface runoff and erosion transported mainly finer textured material, resulting in even coarser texture at the uphill surface substrates.This probably explains the decreasing influence of the substrate proxy variables CORGP and MSAND with time.
Besides such transport effects, the weakening effects of initial surface and soil properties with time reflect the increasing influence of vegetation coverage (see H3).
In summary, our results strongly support H1. Geomorphology (slope degree) and substrate properties like grain size distribution and resulting differences in groundwater levels affect the system strongly and in turn result in fundamental changes in geomorphology all over the area, but with main effects where the slope degree is highest.The developing vegetation cover, however, reduces the influence of these factors in the course of time, leading to an increasing conservation of the structures that developed at the very beginning.Some of our results concerning the influence of surface and substrate characteristics on vegetation structure (H2) match basic assumptions of successional theory such as the decrease of annuals as well as the increase of grasses and especially woody species in the course of time (Cowles, 1899), and the general increase of vegetation height over time (Prach et al., 1997).Despite a scientific history of well over a century (Cutler, 2010) even this is noteworthy.Little effort has been made to integrate and synthesize the large number of case studies in succession to obtain a coherent generalized understanding.This is certainly due to the fact that rarely comprehensive spatial and temporal data on abiotic and biotic aspects of succession are collected.There are studies in the large body of literature which superbly meet a high temporal resolution but with only a very small spatial coverage (e.g.Rebele, 1992).Projects exclusively using chronosequences (probably the most) lack the spatial coverage for a quantitative analysis on vegetation change, both, in the sense of spatial patterns over time and in the sense of temporal pattern in space (Baasch, 2010).Still vast uncertainty exists about process interactions within the plant community (e.g.dispersal, competition) to processes linking the community to its environment (e.g.soil property dynamics) during succession.In view of the above, a particular contribution of H2 and our study is the detailed integration of biotic and abiotic aspects of succession with both temporal and spatial data.
Annual species -most of them being pioneers -were preferentially and expectedly found at places with less organic carbon, which means with low water and nutrient storage capacities (as mentioned before organic carbon might act as a proxy for general substrate properties), and lower groundwater levels.It was found also in a low-pH system that particularly resource-limited places were almost exclusively colonized by pioneers (Baasch, 2010).
In line with basic assumptions, the annual species' portion decreases with time (Fig. 4).However, during the first years of the succession two annual species (Conyza canadensis, Trifolium arvense) were dominant (Zaplata et al., 2011b) and hence their portion increases with increasing vegetation cover (COVTOT) but does not change with COVTOT values beyond 50 % (Fig. 4).The existence of a succession phase governed by Trifolium arvense (Fabaceae) in 2008-2009 means, that many (adjacent) vegetation plots had similar covers of this plant species (Zaplata et al., 2013).The "optimum" shortly before COVTOT = 50% (Fig. 4) might reflect the "preferential cover range", the range of cover which is typical for a given plant species, in this case of T. arvense.Apparently this phase expresses itself in an optimum calendar year effect for the occurrence of Fabaceae (Fig. 7).
In previous studies we show a similar temporal trend for both cover (Schaaf et al., 2013) and regularity of the cover (Zaplata et al., 2013) for the two dominant species.
However, the proportion of annual plants' (to which T. arvense belongs) steady decrease means that the perennials on the whole had disproportionately high growth rates, even during the dominance phases of the two annual species.Besides, the previously mentioned bryophytes' cover increase in 2010/2011 incrementally alters (reduces) cover shares of vascular plant groups.
What has been found for annuals is in many respects also true for herbaceous plants and for Fabaceae (Fig. 7), simply because most of the herbaceous plants at the beginning of the successional series are annuals, and a dominating species of this group (Trifolium arvense) belongs to the Fabaceae family.
For the grasses we found two opposing trends.Our model separates a total cover effect and a temporal effect.The latter indicates a strong tendency to increase the grasses' share on the total community (Fig. 6, Supplement Fig. S6) in contrast to the herbs.Grass species like Agrostis capillaris and Calamagrostis epigejos increased strongly (Zaplata et al., 2011b), and we expect a continued rise in the importance of the grasses.The former reveals that, again in contrast to the herbaceous plants, grasses dominate plots with low total coverage (Fig. 6).Hence, grasses are important contributors on more sparsely covered plots.In such a way they may rule subareas, but not the study system as a whole.This fact supports the assignment of the studied time span to the (first) successional stage of herbaceous vegetation (Zaplata et al., 2013).Only two explanatory variables, total degree of coverage and calendar year, showed a significant relationship in this context (Fig. 6, Table 9, Eq.13).
The same applies for woody plants, with an increasing share on the total community and -unexpectedly -a decreasing probability of finding them on vegetation squares with higher total cover (Table 10).We anticipate a reverse situation in the future.At the latest, when there are several crown layers formed by differently sized trees, their joint covers will prevail over the total covers of pure herb and grass plots.
Based on considerations about resource limitation-driven competition modes by Hara (1993) and Weiner (1990), we hypothesize that our finding indicates that during the observed time span plant growth in the catchment was mostly limited by the soil bound resources, water and nutrients.Under such conditions the tree competition strategy to invest in vertical size does not pay off so much as it does when the vectorially distributed resource light is limiting.As long as this is not the case, trees may often face disadvantages in competition with non-woody plants.
We also found that overall vegetation height tends to be greater on more densely vegetated plots (Table 12).We may interpret that (i) simply as a consequence of particularly poor growth places, (ii) the high proportion of per se low-growing annuals there, or (iii) as an effect of beginning light competition which comes, however, hardly from trees yet.
Our finding that vegetation height tends to be greater when the distance to the groundwater table is higher probably reflects a spurious correlation resulting from the coincidental temporal distribution of dry years co-occurring with ongoing succession.On the other hand there are, however, indications that at least some of this effect might be attributable to the intrinsic characteristics of the system.For example, the first dominant plant species, the aforementioned Conyza canadensis, reached the greatest heights in 2006 (known from orienting measurements on the catchment and an adjacent site), which was the driest year so far in the time series.In the course of time a decrease in size was registered even for some perennial species like Artemisia vulgaris, Hypericum perforatum L. and Hieracium umbellatum, as shown by the presence of higher shoots from the previous growing season.
In conclusion, the results of the separate analyses of plant functional groups (annuals, herbs, Fabaceae, grasses, woody plants) are in line with successional theory, e.g. the decreasing proportions of annuals and herbs and the increasing importance of grasses and woody species.The latter two were under-represented on denser vegetated plots.This seems paradoxical at first sight, considering that woody species are superior competitors.With such results, this study supports a non-neutral view of ecosystem assembly during succession (cf.Nee and Stone, 2003).Our statistical analyses confirm at least in parts hypothesis 2, that surface and substrate characteristics determine the plant species groups to be found initially.Some few substrate characteristics such as the gravel content and the percentage of organic carbon already correlate in different ways with the occurrence of certain plant groups (cf.Baasch, 2010).Expectedly, woody plants show the weakest connection to substrate properties, the weak but significant connection we found between the occurrence of the N-fixing Fabaceae and the soil gravel content (Fig. 7) is plausible as well, as more gravel content often means more unfavourable conditions in terms of water and nutrient supply.On soils with higher gravel content, plants belonging to the Fabaceae family should possess a pronounced competitive advantage: a symbiotic relationship with root nodule bacteria (Rhizobiaceae) allows for the utilization of atmospheric nitrogen improving the supply of this typically limited nutrient.
Introducing vegetation properties into the model for the surface height change (DHEIGHT) clearly confirmed our third hypothesis -the increasing density of vegetation reduces the erosion and therefore the degree of change of terrain surface structure.Remarkably VEGHD, the product of the laser-scanned vegetation height with the laser-scanned vegetation density VEGDENS -a proxy for biomass -performed best as the vegetation-related explanatory variable.The effects of all other substrate and geomorphology variables stayed virtually the same as in the DHEIGHT model without any vegetation variables used in the context of our first hypothesis.As Fig. 9 shows, the terrain-stabilizing effect of vegetation increases in the course of time, even for the same values of VEGHD.This may be interpreted as the result of positive (reinforcing) feedbacks between vegetation and substrate (vegetation establishment stabilizes the surface which supports further vegetation establishment) while vegetation growth itself appears to be also driven by a positive feedback as can be seen from the exponentially shaped calendar year effect on vegetation height in

Conclusions
The initial geomorphology in terms of position along the slope and the substrate properties, organic carbon and medium sand content have, together with groundwater levels, an effect on the surface and its formation as expressed by local rill density, relief energy, and surface height change.Expectedly, the influence of initial substrate properties diminishes with time (H1).The occurrence of plant functional groups (annual, herbaceous, grass-like, woody, Fabaceae) is connected to initial substrate properties; however this connection does not seem to be particularly close (H2).Significant relationships could be identified only with organic carbon content and, in the case of Fabaceae, with gravel content.Groundwater level and the local degree of total plant cover more often turned out as important explanatory variables for the appearance of a given plant functional group as well as for local mean vegetation height.For vegetation density no significant variables were found.In general, the time-dependencies identified for the occurrence of different plant functional groups are in line with widely accepted views on succession.The product of vegetation height and density, a proxy for vegetation biomass, showed an increasing counteracting effect on surface height decrease (erosion), thus confirming H3.
Our analysis reveals dynamic processes driving the ecosystem's development.We find positive (reinforcing) feedback mechanisms like the tendency of rills to deepen.However such pure substrate or geomorphology induced effects weaken with time, while vegetation establishes increasingly and differentiatedly introducing new (currently positive) feedbacks.These consist in mutually reinforcing stabilization of surface and vegetation, but also in vegetation growth processes.The analysed period appears to represent the transition from a geo-hydro-system towards a bio-geohydro-system as defined in Elmer et al. (2013).We strongly believe that continued and broad research in that particular catchment will substantially contribute to an improved understanding of ecosystem functionality -its drivers, mediators, regulators, preservers, and antagonists.

Fig. 1 .
Figure 1 1092 Fig. 1.Map of the Chicken Creek area (with Gauss-Krüger zone 5 coordinates) showing locations of the data sources (vegetation plots, laser scanner positions, groundwater gauge positions) relevant for this study.Vegetation plots are identified by "row" (A-T) and "column" (1-7).In this work only rows A-Q were used.The variable DISTA as shown in the figure indicates the distance of a point from row A. Ground heights a.s.l. are given as mean values from the laser scan in September 2008 around the centre points of vegetation plots A6, F6, L6 and Q6.The height profile along row 4 is given (top right) in the map.It shows the ground elevation as obtained from the laser scans in September 2008.

Table 1 .
Names, data source and description of the variables used in this study.

Table 2 .
Summary statistics of the variables used in this study (see Table1for variable explanations).DHEIGHT, indicating a change between subsequent years, is given always at the start year.The soil-related variables CORGP, GRAVCONT, and MSAND were only measured once in 2005.DISTA does not change with time and is presented for the time when the catchment's establishment was completed.The variable COVTOT is shown in Table3together with other vegetation properties.

Table 3 .
Minimum, mean and maximum covers in % per plant species group and year of survey.The allocation of the single plant species to the groups is shown in Table14.COVTOT is the total degree of coverage by plants (see Table1).This variable has larger values than the cover sum of annual and perennial, or grass-like, herbaceous, and woody plants, because it also includes the bryophytes.