Uncertainty in parameterisation and model structure affect simulation results in coupled ecohydrological models

In this paper we develop and apply a conceptual ecohydrological model to investigate the effects of model structure and parameter uncertainty on the simulation of vegetation structure and hydrological dynamics. The model is applied for a typical water limited riparian ecosystem along an ephemeral river: the middle section of the Kuiseb River in Namibia. We modelled this system by coupling an ecological model with a conceptual hydrological model. The hydrological model is storage based with stochastical forcing from the flood. The ecosystem is modelled with a population model, and represents three dominating riparian plant populations. In appreciation of uncertainty about population dynamics, we applied three model versions with increasing complexity. Population parameters were found by Latin hypercube sampling of the parameter space and with the constraint that three species should coexist as observed. Two of the three models were able to reproduce the observed coexistence. However, both models relied on different coexistence mechanisms, and reacted differently to change of long term memory in the flood forcing. The coexistence requirement strongly constrained the parameter space for both successful models. Only very few parameter sets (0.5% of 150 000 samples) allowed for coexistence in a representative number of repeated simulations (at least 10 out of 100) and the success of the coexistence mechanism was controlled by the combination of population parameters. The ensemble statistics of average values of hydrologic variables like transpiration and depth to ground water were similar for both models, sugCorrespondence to: S. Arnold (sven.arnold@ufz.de) gesting that they were mainly controlled by the applied hydrological model. The ensemble statistics of the fluctuations of depth to groundwater and transpiration, however, differed significantly, suggesting that they were controlled by the applied ecological model and coexistence mechanisms. Our study emphasizes that uncertainty about ecosystem structure and intra-specific interactions influence the prediction of the hydrosystem.


Introduction
In semiarid environments water is not only a scarce resource, water availability also varies greatly in timing and magnitude.Both natural ecosystems and people have to adapt to these conditions, and often they share the same water source.Thus, water management of the water source might influence natural ecosystems, but also inversely, the management of vegetation might affect the water fluxes.In order to understand, what implications human development in semiarid regions has, models are required that help investigating the effect of management actions.Such models need appropriate description of both ecological and hydrological processes.
A great deal of work in ecohydrology has already been dedicated to understanding mechanisms, by which a variation in water availability influences vegetation patterns.Much of this work is based on considering single plant species, and comparing expected water stress-levels in different environments.Therefore, these models cannot consider inter-specific competition or coexistence.However, research dealing with biodiversity and species-co-existence suggests that particularly fluctuations of environmental signals might favour co-existence (D'Odorico et al., 2008).Hence certain levels of variance of water availability could also be a driver for maintaining multispecies plant communities.Moreover, diverse ecosystems are thought to be more resilient to disturbance and should thus react differently to extreme conditions than single species ecosystems.Hence coexistence mechanisms might be important ecosystem processes shaping plant-water interactions in water limited environments, which motivates the need for multispecies ecohydrological models.Such a model is developed and applied in this paper.
Ecological modelling has different approaches to describe multi-species plant communities.One way is spatially explicit individual-based modelling (DeAngelis and Gross, 1992;Grimm and Railsback, 2005), representing a bottomup approach.Here, plant communities are described as systems of interacting plant individuals responding to their environment.This approach is particularly powerful when specific systems are to be analysed.The respective models, however, are often complex that makes parameterisation a challenge (lots of parameters) and hampers generalization (adjustment to a specific case vs. principle understanding, transferability).To gain principle understanding of the interplay between water resources and vegetation and the response of environmental variability along ephemeral rivers is central for the present study.Therefore, we follow a topdown approach, i.e. we use a multi-species population dynamical model (Kot, 2001) to describe the plant community in an aggregated way but explicitly consider the species' competition for water.The population dynamical parameters summarize all relevant effects caused by processes at the individual scale (e.g.plant growth and mortality, response to disturbances, type and strength of competition, seed dispersal) (Fahse et al., 1998;Frank and Wissel, 2002;Heinz et al., 2005;Moorcroft, 2003;Ovaskainen and Hanski, 2002;Ovaskainen and Hanski, 2004).Most ecohydrological models work at the population scale (Camporeale and Ridolfi, 2006;Porporato et al., 2001;Ridolfi et al., 2000;Rodriguez-Iturbe et al., 2001).Direct parameterisation of population models is sometimes impossible as this requires long-term observation of species abundance, which is not always available.
Generally, both population and hydrological models can be developed with varying levels of complexity.In order to keep a coupled model manageable, the level of model complexity needs to be appropriate regarding the desired predicted variable but also regarding the available data.And there has to be a strategy how the model should be parameterised.
In this study, we address this parameterisation problem by using pattern-oriented model calibration, in that we adjust species parameters such that the resulting model reproduces the observed coexistence.Models have been parameterised based on information of presence or absence of plant species before (Laio et al., 2001;Rodriguez-Iturbe et al., 1999).When the existence criterion will be extended to sev-eral species it is called coexistence, and also observed coexistence has been used to evaluate (at least qualitatively) the validity of ecohydrological models.In doing so, researchers put their models to a strict test, since modelling coexistence is comparatively difficult (Arora and Boer, 2006;Clark et al., 2007).A given model only allows for coexistence, if its structure and parameters meet strict conditions, which provide for the required relation of trade-offs.A number of mechanisms can be invoked fostering coexistence in models, such as ecological niches (in time and space) and tradeoffs (Chesson, 2000;Clark et al., 2007).Ecological theory also indicates that the variability of an environmental signal, such as resources or disturbance regimes, influences biodiversity.According to the Intermediate Disturbance Hypothesis (Connell, 1978;Huston, 1979), moderate levels of environmental fluctuations can enhance both biodiversity and resilience (D'Odorico et al., 2008).So far, such studies have dealt with uncorrelated, random environmental signals.Examples are given for random water table (Ridolfi et al., 2007) and climate fluctuations (Rodriguez-Iturbe et al., 1999), or environmental disturbances induced by fires (Higgins et al., 2000;van Wijk and Rodriguez-Iturbe, 2002).However, many hydrologic time series are characterized by auto-correlated and longterm-memory processes (Hurst, 1951;Montanari et al., 1997), particularly in arid environments.This directly leads to the question of the role of this autocorrelation, that is the duration of a disturbance event (water stress, disruptive flood), for the functioning of the ecohydrological system.Moreover, studies usually consider only one consequence of an environmental signal.However, the same signal, for example rain, may interact with the system in multiple ways.A strong rain event might recharge the water storage for plants, but at the same time, the storm might destroy part of the vegetation.Thus the event acts on both, mortality and growth, but possibly not in the same fashion.Such combined effects are not fully understood so far.In this work we wish to investigate both of these issues, based on the example of an ephemeral river in Namibia.This allows for testing the adequateness of the Intermediate Disturbance Hypothesis in the context of ecohydrological systems along ephemeral rivers.
The middle section of the ephemeral Kuiseb River in Namibia is a representative example of an environmental system with ecohydrological feedbacks and need for management.Previous studies indicate that the development of riparian vegetation depends on the subsurface water storage (alluvial aquifer) which is recharged by intermittent floods.At the same time, strong floods lead to uprooting of riparian vegetation and increased mortality.There is negligible rainfall in this part of the river, the floods originate in the upper reach, and depend both on the rainfall regime and small scale farm dams.In this study, we aim to build a model that allows understanding, how the flood regime interacts with the riparian ecosystem and the resulting transpiration loss and aquifer storage.Little data is available regarding the ecosystem.We therefore rely on conceptual models both Hydrol.Earth Syst.Sci., 13,[1789][1790][1791][1792][1793][1794][1795][1796][1797][1798][1799][1800][1801][1802][1803][1804][1805][1806][1807]2009 www.hydrol-earth-syst-sci.net/13/1789/2009/ for ecosystem and aquifer.In order to address structural uncertainty, we select three models, with increasing degree of complexity of the ecological model.We attempt to parameterize these models based on the scarce available information, namely the fact that three species coexist and some knowledge about their maximum transpiration rates and rooting behaviour.Our investigation shows that different coexistence supporting mechanisms can be invoked, depending on the assumed conceptual model.While the distribution of mean hydrologic variables (groundwater level and transpiration) was similar in all models, their variability depended both on the model structure and the parameters sets.This points at the difficulty to parameterise an ecohydrological model in real world applications.However, our model gives clear indications, what measurements are most effective for improving the necessary process understanding.
2 Methods and materials

Study site
The study site covers an area of approximately 18 km 2 and is located in the Kuiseb catchment (∼15 500 km 2 - Jacobson et al., 1995) in Namibia (Fig. 1).The Kuiseb River arises from the Khomas Hochland (∼2000 m in elevation) and runs westward through the escarpment into the Atlantic Ocean.The rainy season is during the Southern Hemisphere summer between January and April (Henschel et al., 2005).Most of the rain falls in the upper reach of the catchment (Khomas Hochland).This study is concerned with the arid middle reach of the Kuiseb River, where rain is exceptional, and water arrives mainly during the floods in the ephemeral river channel.Near this channel, riparian vegetation has established.Although the channel does not contain water for most of the year, it supplies a shallow aquifer with water during times of flood and thus creates a living environment for riparian vegetation.The flood is influenced by upstream farm damns and the ground water table is influenced both by plants and human consumption.

Ecosystem
Vegetation around the river channel consists of 80% of only three coexisting species: Camel Thorn (Acacia erioloba), Ana Tree (Faidherbia albida) and Wild Tamarix (Tamarix usneoides) (Theron et al. 1980).All of them depend on the infiltration of flood water, with slight differences in strategies.Schachtschneider and February (2007) investigated the water use strategies of all three species by using isotope methods.They found that both Camel Thorn and Ana Tree use a mixture of ground-and soil water, and Wild Tamarix uses water from the unsaturated zone, originating from flood and also fog water.The known differences between the three species are in their phenology (time of leaf shedding), maximum transpiration and growth rates (see Table 1).Besides supplying vegetation with water, floods in the river channel have also a destructive component.Small trees are usually washed out by strong floods.The latter makes slow growing trees vulnerable for large floods for longer time.

Hydrosystem
The study site is located in a hyperarid area with mean annual rainfall less than 20 mm and mean potential evaporation of 1700 to 2500 mm (Botes et al., 2003).The shallow alluvial aquifer consists of sand and is embedded into impermeable granite (Dahan et al., 2008;Morin et al., 2009;Schmidt and Plöthner, 1999) (Fig. 2).Its thickness and width vary along the river.The alluvial aquifer is recharged by temporary floods that are caused by rainfall in the upper Kuiseb catchment (Khomas Hochland).Volume and duration of the resulting floods vary strongly (Fig. 3).Larger floods burst over the limits of channel bed, leading to inundation of the river banks.At the same time, about 90% of the floods run dry within the Kuiseb middle section under study here.This shows the comparatively large role of infiltration.The dynamics of flood water infiltration were investigated by Dahan   (Canadell et al., 1996), c (Dalpe et al., 2000), d (Stave et al., 2005), e (Moser, 2006), f (Schachtschneider and February, 2007), g (Coates Palgrave, 1983), e (Timberlake et al., 1999), h (Wickens et al., 1995).et al. (2008).Their studies show that, during a flood, the water content of the unsaturated layer only increases up to the twofold value of the field capacity and that the infiltration rates across different flood events are very similar.Further, above a certain flood stage threshold, it is the flow duration and not the flood height that controls the recharge amounts.

Hydrological Model
We modelled the hydrological processes along an ephemeral river with shallow aquifer.Figure 2 gives a sketch of the hydrological unit modelled.We modelled a representative river-valley segment of 60 km length and a constant width of 300 m.Hence we considered total fluxes over the entire surface area of the segment, which is A seg =18 km 2 .The water balance for this segment is written as where WS(t) is the sum of change in unsaturated ( S unsat (t)) and ground water storage ( S GW (t)).The storage in the unsaturated and ground water layer was calculated as •φ for the ground water storage, (2b) where θ (t) is the water content (m 3 /m 3 ) of the unsaturated zone, which ranges between 0 and porosity φ (Table 2) and h GW (t) is the ground water level.The depth to ground water z unsat (t) is where h WS =15 m is the total depth of the alluvium.In our simulations we fixed the initial value of ground water depth to z unsat (t=1)=5 m.
Hydrol   4).We used the Hydraulic Properties Calculator of Saxton and Rawls (2006) to estimate the volumetric water content at permanent wilting point (θ PWP ), field capacity (θ FC ) and the porosity (φ).For this study we assumed the soil texture class of the alluvial fill to be sand with an average grain size distribution of 8% gravel, 90% sand, and 2% clay.

Soil
Flood generator Value (FARIMA) The change in unsaturated storage was calculated as with I (t) denoting the infiltration, GWR(t) the ground water recharge, and T unsat (t) the transpiration from the unsaturated storage.The infiltration to unsaturated soil is based on the results of Dahan et al. (2008) who concluded that infiltration fluxes are limited by a flux-regulating mechanism at the top of the unsaturated zone, independent of the flood height.They suggest a time constant infiltration rate of Q I (t)=1 cm/h, which is 2400 m 3 /d ha.Therefore, the infiltration depends only on flood duration D(t) (Eq.22) and the specific infiltration flux Q I (t): Flood duration is calculated as a function of flood volume in Eq. ( 22) (see Sect. 2.3).The ground water recharge depends on the water content θ(t) of the unsaturated layer: where S FC (t) is the water volume in the unsaturated zone corresponding to the water content at field capacity (θ FC (t)=0.061).The transpiration is composed of transpiration from unsaturated layer and ground water.The transpiration from the unsaturated layer is the sum of the transpiration from individual species T unsat,i (t): where S PWP (t) is the water volume in the unsaturated zone corresponding to the water content at permanent wilting point (θ PWP (t)=0.015).For plants where the roots reach the groundwater, transpiration originates from both the unsaturated and the saturated zone.The unsaturated part is calculated as where T WS,i (t) is the transpirational demand for each species (Eq.11), V unsat,i (t) the water volume in the unsaturated storage and V WS,i (t)is the total water volume (unsaturated and ground water) that can be reached by plant roots of species i: where V GW,i (t) is the ground water volume available to plant roots.The water in the unsaturated storage available for transpiration of species i depends on its rooting depth z r,i (t): Note, that for the purpose of keeping the model simple we neglected any age structure in the ecological model (see Sect. 2.4).Consequently, the rooting depth does not depend on the age of a (sub)population.The transpirational demand for each species (T WS,i (t)) is a linear function of the green biomass G i (t) (see Sect. 2.4) with an upper boundary given by the potential evapotranspiration (PET): The PET was estimated using the Penman-Monteith Equation for both the flooding and the dry season.The transpiration per green biomass Q T ,i (t) of each species is derived from measurements of Bate and Walker (1991) and is summarized in Table 3.The change in ground water was calculated as where Q GW (t) is the ground water flow and T GW (t)the transpiration of all species from ground water (Eq.15).The ground water flow is where Q In is the ground water inflow from upstream, Q Out (t) the ground water outflow downstream, Q L the lateral ground water inflow, and Q V the vertical ground water outflow to the bedrock.Q In , Q L and Q V are assumed to be constant over time (Table 7).Q Out (t) was calculated by Darcy's Law, as: with k f denoting the hydraulic permeability of the ground water layer, h(t) the hydraulic gradient between the inlet and outlet of the modelled aquifer segment, and A GW the cross-sectional area of the ground water layer.The transpiration of all species from ground water is the sum of individual species transpirations T GW,i (t): where V GW,i (t) is the ground water that can be reached by plant roots of species i.
In the water balance described above, we neglected two processes: precipitation and evaporation.The first is very low at the study site (23.8 mm/year at Gobabeb Research Centre -Schulze 1969).The second is only active during flooding, which is only a few days per year.The effective depth of direct evaporation from bare soils was assumed to be 1.5 m and can be considered as non active soil layer above the alluvium.

The stochastic flood generator
The flood volume V Flood (t) was generated by a fractional autoregressive moving average (FARIMA(p, d, q), p, q∈N ) model with symmetric α-stable (SαS, α∈(1,2)) innovations (Kokoszka and Taqqu, 1995;Stoev and Taqqu, 2004).The FARIMA(p, d, q) model generates time series with both short-and long-term dependence structures that are present in many hydrologic processes (Hurst, 1951;Montanari et al., 1997).We used the algorithm presented in (Stoev and Taqqu, 2004) to generate time series with given short-and long-term memory.The short term dependence structure is determined by the real polynomials X p and q of degree p and q.The autoregressive part of FARIMA is represented by the coefficients of X p , where X 1 (λ)=1−0.192λand λ is a random number drawn from a normal distribution with mean 0 and standard deviation 1.The moving average part is represented by the coefficients of q : with 1 (λ)=1−0.8969λ.The long term behaviour is governed by d that is an arbitrary fractional real number: The relationship between d and the Hurst-Exponent H is as follows: The value of H varies between 0 and 1, an H of 0.5 means absence of long term memory or white noise.Values lower than 0.5 correspond to negative dependence; however, these are rarely encountered in the analysis of hydrologic data (Montanari et al., 1997).Typical values of H range between 0.7 and 0.8 (Hurst, 1951).Hence, for our study, we assumed H to be 0.75 (with α=1.99 and d=0.25), and p=q=1.
The time series were generated with FARIMA(p=1, d=0.25, q=1) and adjusted to the observed mean annual flood volume µ Flood =3 269 000 m 3 , and thus yielding Flood duration was found to be related to flood volume.Therefore we performed a linear regression between the measured flood duration and the corresponding logarithmic flood volumes from 1981 to 2006.The derived best fit (r 2 =0.9) was given by and used in the following to calculate the flood duration.

Ecological model
The ecological model aims to describe the dynamics of the plant community consisting of the three tree species of interest in the river-basin of the Kuiseb in relation to the availability of water as jointly utilized resource.Each tree species is characterized by its biomass in the river-valley segment.In order to address important processes of the plant community dynamics and their response to the hydrological system in an adequate way, biomass of a species is differentiated into green (G) and reserve biomass (R) similarly as (Muller et al., 2007), who termed R after (Noy-Meir, 1982).The green biomass describes all the parts of a plant, which perform photosynthesis, while the reserve biomass covers all parts of the plant that are not photosynthetically active, like woody parts and roots.The dynamic of G is driven by seasonality (phenology) and short-term water stress.The process of photosynthesis performed by G depends on the availability of water (transpiration, see Sect.2.2) and results in the production of organic carbon, which maintains both green and reserve biomass.The dynamic of R occurs on a longer timescale and reflects the long-term history of the ecohydrological system.The model is applied at a seasonal time scale, thus dividing the year in two halves: the season when floods occur (southern hemisphere summer) and the dry season.During the seasons, when the plants are photosynthetically active, the green biomass G i is modelled as where G i (t) and G i (t−1) are the green biomass in this and the previous time step of species i, with units of t/ha, R i (t−1) is the reserve biomass in the previous time step, w G,i (t) is the conversion rate from reserve into green biomass (Eq.26), and ε i (t) is the unitless water stress function (Eq.27), ranging from 0 for no water stress to 1 for complete water stress.The latter two terms, w G,i (t) and ε i (t), are functions of the available amount of water (Eqs. 26,27).The first term of Eq. ( 23), (1−ε i (t))•G i (t−1), denotes the leaf shed due to water stress, while the second part, w G,i (t)•R i (t−1), denotes the growing of leaves on the existing reserve biomass, assuming that the required Carbon of the reserve biomass was already accumulated in the buds during the previous season.
Depending on the complexity of the model, we either assume no phenological differences between the species (model A), or we include the known differences in phenology.In the first case, Eq. ( 23) applies to all species at all times.In the latter case, some species are dormant during a particular season (model B and C). Green biomass during the dormant season was calculated as: where ls i is the unitless leaf shedding factor and ranges from 0 to 1 of species i. ls i =0 corresponds to no leaf shedding at all, and 1 to complete leaf shed.Usually, leaf shed is not complete, so ls i takes a value between 0 and 1.The formation of reserve biomass takes place at the end of each season t: where f r i (t) is the unitless flood resistance of species i and ranges from 0 to 1 (Eq.28).It denotes the vulnerability of a given species to being uprooted and washed away by a flood of given magnitude.f r i (t)=0 corresponds to complete removal of reserve biomass by the flood.In the dry season, f r i (t) is set to 1.The parameter m R,i denotes the mortality of the reserve biomass, and w R,i the growth rate of reserve biomass.Both are constant over time and unitless.The first part of Eq. ( 25), [1−m R,i (1+ε i (t))]•R i (t−1), denotes the amount of reserve biomass remaining after mortality and response to water stress.Note that the total mortality increases when ε i (t)>0.The second part, w R,i •G i (t), corresponds to growth of reserve biomass, based on the photosynthesis performed by the green biomass G i (t).In our simulations we fixed the initial values of green and reserve biomass to G i (t=1)=0 t/ha and R i (t=1)=0.1 t/ha.
In our model, favourable periods of growth in the green biomass G can markedly increase the reserve biomass R, whereas unfavourable periods reduce G fast, but R only slowly.In his paper about the multispecies competition in variable environments, Chesson (1994) called this the storage effect, which "is a metaphor for the potential for periods of strong positive growth that cannot be cancelled by negative growth at other times".The storage effect is enhanced by the parameter w R,i (Eq.25).
The three parameters conversion rate w G,i (t), water stress ε i (t), and flood resistance f r i (t) are characteristics of the tree species that are dynamically linked to the hydrosystem.The conversion rate from reserve to green biomass, w G,i (t), is described by a sigmoid function that depends on the water volume in the alluvium that can be reached by the plant roots (V W S,i (t)) (see Sect. 2.2, Eq. 9) and the total reserve biomass of the ecosystem in the previous time step (R total (t−1)= where a i , b i and c i are the shape parameters of the sigmoid function, and depend on species i.The dependence of w G,i (t) on accessible water volume V WS,i (t) and total reserve biomass R total reflects the intra-and interspecific competition between the three plant species for water, although in an aggregated and non-spatial way.The water stress function ε i (t) was calculated as where V stress,i is the water volume in the alluvium reachable by plant roots that leads to water stress in the population of species i and V PWP,i is the water volume within the reach of plant roots that is no more extractable by plants.It is speciesspecific because it depends on the species root depth Eqs.(9, 10a).V Stress,i is also a species-specific parameter: the lower V Stress,i the more drought tolerant is this species.The flood resistance f r i (t), describes the capacity of the vegetation to withstand a flood without being uprooted and washed away.It reduces the reserve biomass, which is assumed to be built at the end of season, and only applies during the flood season (Eq.25).We modelled it as a linear function of the flood volume (V Flood with unit m 3 /ha), which was generated by Eq. ( 21).f r i (t where f i and g i are species specific shape parameters.When the flood volume is below V low,i the flood resistance is 1, the flood is minor and the species population does not suffer additional mortality induced by flood.Above the flood volume of V high,i the flood resistance is 0, i.e. the species population is completely washed away.

Model versions
One aim of this study was the analysis of model complexity with regard to model output.Therefore we investigated several model types that differ in complexity regarding the representation of the ecosystem.Since little is known about the ecological parameters in the Kuiseb River, any model would be comparatively simple.We compared three model version of the same area.Table 4 gives an overview about the model differences.
In the first model (A) we neglected phenology, all species were evergreen, but differed other traits like maximum transpiration rate.Generally, leaf shed is only partial for all species, thus it suggests itself to neglect seasonal variation.In the second version B, we included the observed species specific phenology of Camel Thorn and Ana Tree (Table 1).For this we added two parameters (ls Cam and ls Ana ), which increased the degree of complexity (model type B).Finally, in model C, we included more knowledge regarding the dif- ference in flood resistance between species, thus allowing the parameter f r i to be species specific.In summary, model type C included the most ecological information, strongest constraints and a mortality that is not only stochastic but also depends on the hydrosystem.In each model application we compared, if the model was able to reproduce the observed coexistence of three species.To achieve this goal we parameterised the models accordingly, as pointed out in the next section.

Parameter sampling
Depending on the model version, the ecological model contained 23-29 parameters.Table 5 gives an overview of those parameters together with their physical range.We used Latin hypercube sampling in order to identify parameter sets, which lead to the observed coexistence of three species.This was performed for each model version separately.Only the ecological parameters were calibrated, the hydrological parameters were fixed to the values indicated in Table 2.
We constrained the parameter space qualitatively according to the available ecological information summarized in Table 1: the root depth was largest for Camel Thorn, followed by Ana Tree and Wild Tamarix.Further, we assumed that the growth rate of reserve biomass can be derived from wood density, that is, the larger the wood density the smaller is w R,i .Hence, reserve biomass growth rate was largest for Ana Tree, followed by Wild Tamarix and Camel Thorn.Additionally, we checked the sampled parameter sets for plausibility: for the shape parameters a i , b i , and c i we allowed only combinations that lead to w G,i (t)=0, if R total (t−1) ≤0 in Eq. ( 26).
The sampling procedure is illustrated in Fig. 4. For each sampled parameter set ( ι ) we run the model 100 times.We than counted the number of runs, where all three species coexisted and defined the probability of coexistence (P 3,ι ) for the parameter set ι as follows: where # B(n=3) is the number of flood realisations that led to coexistence of all three species.P 3,ι gives an indication how robust the modelled coexistence was.If P 3,ι is small, the Hydrol.Earth Syst.Sci., 13, 1789-1807, 2009 www.hydrol-earth-syst-sci.net/13/1789/2009/ parameter set ι only led to coexistence under very specific flood conditions, while a P 3,ι near 1 indicates that the parameter set led to coexistence in almost all flood realisations with the same stochastic properties.
We defined coexistence based on the following criterion: the average reserve biomass during the last 1000 years must exceed the reserve biomass necessary to maintain 10 adult individuals of average size of each species.The method for deriving the number of individuals is described in the Appendix A.

Analysis of the ensemble models
For analysis of the model results, we used ensemble statistics of hydrological variables of interest and each parameter set ι (with P 3,ι as indicated).The statistics were only performed on the last 1000 years (2000 time steps) of each simulation, in order to avoid the influence of initial conditions.
The expected ensemble mean for parameter set ι of the variable of interest (for example total ecosystem transpiration) was calculated as follows.We first calculated the time series means of the variable of interest, for each simulation that led to coexistence with the same parameter set ι .Secondly, we calculated the ensemble mean of the obtained set of time averages.
We only calculated the time average for the subset of η 3 simulations, which led to three species coexistence.The statistics was performed on at least η 3 =10 simulations.If necessary, additional forward simulations were run in order to obtain 10 simulations with coexistence.Each time average of the variable of interest ( Vι,η ) is calculated as where V t is the value of the hydrological variable of interest at time step t, and η the number of the model realisation.
This led to a set of η 3 time averages for the variable of interest.Based on this set we calculated the ensemble mean of η 3 realisations, which is We proceeded similarly, to obtain the ensemble average of the coefficient of variation of the hydrologic variable.We calculated the dimensionless coefficient of variation (CV V ι ,η ) for the time series: where σ V ι,η denotes the standard deviation within the time series of the variable of interest.Based on this we calculated the ensemble mean of η 3 realisations, which is

Forward simulations with changed flood regime
After finding ensembles of suitable parameter sets, we tested how models behaved for changed flood conditions.For this we selected those parameter sets which led to coexistence, and run them again with changed flood regime.We changed the long term memory of the flood generation algorithm, by decreasing and increasing the Hurst exponent (Eq.20).We grouped the forward simulations into those performed with parameter sets of weak robustness (0.1≤P 3,ι ≤0.5) and elevated robustness (P 3,ι >0.5).

Results
Table 6 shows in per cent how many of the 150 000 sampled parameter sets led to P 3,ι ≥0.1 for models A, B and C. In all cases, the number of parameter sets that allowed for coexistence of all species is very small (less than half percent in all cases).Furthermore, coexistence was modelled for more parameter sets in models A and C, compared to B: the total number of parameter sets leading to P 3,ι ≥0.1 for model A and C was about 20 times (both around 0.2%) larger than model B (0.009%, Table 6).Model B was not subject to further investigations because there were no parameter sets leading to elevated robustness of three species coexistence with P 3,ι >0.5.
In Fig. 5 we plotted histograms of the achieved probabilities of coexistence (P 3,ι ≥0.1) for model A and C.These histograms give an impression how robust the modelled coexistence was for the different models.The skewness γ of both histograms indicates that most parameter sets showed little robustness (γ A =1.5 and γ C =1.7).Also, for model A the number of robust parameter sets was larger.For example, consider only parameter sets with 0.1≤P 3,ι ≤1: in model A 14.3% of those had P 3,ι >0.5, but in model C only 4.2%.
In order to show how models A and C differ hydrologically we compared the distributions of the ensemble means of hydrologic variables for parameter sets ( ι ) with P 3,ι ≥0.1.In Fig. 6 we plotted histograms of the ensemble average of total transpirations (left) and depths to ground water (right).In model A the transpiration was larger (median 161 mm/year) than in model C (median 148 mm/year).In contrast, the depth to ground water was similar for both models (median A: 7.38 m, C: 7.63 m).The difference for model A and C becomes apparent when comparing the extremes of depth to ground water.In model A the ground water was more often modelled close to the surface (0.25 percentile was 5.93 m) than in model C (0.25 percentile was 6.63 m).The opposite is true for deep ground water tables (0.75 percentile in model A was 11.80 m versus 9.74 m in model C).
In Fig. 7 we plotted histograms of the ensemble means of CV for total transpiration and depth to ground water.While models A and C differed little with regard to the distributions of the ensemble averages of transpiration and depth to ground water, they were much different with regard to the distributions of the time fluctuations of these variables.In model A the time fluctuation in transpiration was much lower (median 0.258) than in model C (0.799).Less pronounced was the difference in the variation of ground water depth, which was also smaller in model A (median 0.025) than in model C (median 0.084).
Next, we investigated, if increase in robustness was related to similar parameter sets and similar hydrological conditions.In other words, are all robust parameter sets just small variations of a similar model, or are they completely different?For this, we looked at both the modelled hydrology and the difference between parameters.In Fig. 8 we plotted the medians of transpirations and ground water depth corresponding to the probabilities of coexistence (P 3,ι ).In Fig. 9 we plotted the medians of CV of transpiration and ground water depth corresponding to the probabilities of coexistence (P 3,ι ).Both, Figs. 8 and 9, suggest that in model C a weak relationship existed between the robustness of the parameter sets (P 3,ι ) and transpiration (r 2 =0.34) and ground water table (r 2 =−0.38).Also, a weak relationship existed between P 3,ι and the CV of transpiration (r 2 =−0.27) and ground water table (r 2 =0.14).No such relation existed for model A. Figure 10 gives an impression how robustness of the parameter sets was related to the similarity of four parameters in model C: the root depth (z r,i ), the growth rate of reserve biomass (w R,i ), the mortality of reserve biomass (m R,i ), and the shape parameter c i of the conversion rate from reserve to green biomass.The plots show that no relationship between robustness of the parameter sets and parameter similarity existed.The same holds for model A.
In Figs.11 and 12 we plotted typical time series of the reserve and green biomass, the flood volume and the depth to ground water.These time series allow insight into the driving coexistence mechanisms in model A and C. In model C the biomass and ground water was more affected by the flood (Fig. 12a) than in model A (Fig. 11a).In model C, two alternating states existed.One state was associated with high prevalence of Camel Thorn and Wild Tamarix, small floods and deep ground water table (e.g.year 600-750 in Fig. 12).The other state was associated with high prevalence of Ana Tree and Wild Tamarix, strong floods and shallow ground water table (e.g.year 850-1000 in Fig. 12).In all parameter sets of model C Ana Tree was characterized by a  larger vulnerability to flood disturbance than Camel Thorn and Wild Tamarix (Fig. 12b).Model A showed different dynamics.In model A the green biomass and the ground water remained constant after initial fluctuations (Fig. 11).The time series of each species reserve biomass were synchronized with small and frequent disturbances by the flood.
Figure 13 shows how model A and C were affected by a changed long term memory of the flood volume.The rel-ative frequencies refer to the previously identified parameter sets with low robustness (0.1≤P 3,ι ≤0.5, Fig. 13a) and elevated robustness (P 3,ι >0.5, Fig. 13b).In model C decrease of long term memory decreased species coexistence.This effect was even stronger for the robust parameter sets (Fig. 13b).In model A three species coexistence was little affected by change of the long term memory of the flood, and independent of the robustness of the parameter sets.

Discussion
We applied three ecohydrological models that differ in the amount of included information, and structure.Differences particularly concerned the functional response of the plant species to the hydrosystem along the ephemeral Kuiseb River.We assessed these models regarding their ability to predict coexistence of the three species as was observed in reality.This strategy of pattern-oriented modelling (see e.g.Grimm et al., 2005 and references therein) has been used to model coexistence before.In our study, only two of the three models allow for robust coexistence of all three species.Further, in both models only few parameter sets reproduce coexistence.This is in line with the classical competition models from ecology (e.g.Lotka, 1925;Volterra, 1926).These models also reveal that species coexistence only emerges if certain restrictive conditions are met by the model parameters.As a result, the parameter combinations found to be appropriate are sparse given the entire parameter space.The comparison between observed and simulated patterns acts as a filter, which allows us to identify, whether a given model structure and parameter combination allows coexis-tence.In this study, only models A and C allow for robust coexistence.They describe two different coexistence mechanisms for different levels of detail.In model A, species are found to co-exist only, if they have access to different water storages, depending on their root depths (Fig. 11b).Camel Thorn has access to deep ground water and does not compete with any other species.On the other hand, the roots of Ana Tree and Wild Tamarix can only reach the unsaturated layer.Hence, only these two species compete for water in the unsaturated layer.Their coexistence is driven by the trade-off between growth rate of reserve biomass (w Ri ) and water stress (ε i ), both influencing green biomass and, hence, transpiration demand of the individual species (see Eq. 11).Ana Tree, for instance, has the larger growth rate, but is less water stress resistant.Therefore, coexistence in model A is based on both niche partitioning and trade-offs.
In model B this sensible balance is broken, by introducing the (observed) phenology.The phenology of Ana Tree in model B reduces the growth period to one season whereas the direct competitor, Wild Tamarix, is evergreen and uses the water resource all year.This provides Wild Tamarix with an advantage in the competition over Ana Tree.In other words, inter-specific competition is enhanced in Model B with the effect that coexistence of all three species is not possible anymore.This is in accordance with the classical competition theory (see above).Note that this also indicates that integrating more knowledge in a model does not automatically lead to more realistic modelling results.On the other hand, models can give satisfactory results, but maybe for the wrong reason.Effects may be neglected which can play an important role under different management or climatic conditions.
In model C, another coexistence mechanism is enabled, only by allowing for species specific vulnerability to the flood.Thus, as opposed to models A and B, the flood has differential influence both as a water resource and via the destructive impact of the flood; the latter acts directly as an environmental disturbance on the plant species and favours flood resistant species during periods of strong floods.This can compensate the disadvantage of being less competitive than other species in other respects and, hence, can mediate Fig. 10.Parameter space of root depths (z r,i ), growth rates of reserve biomass (w, R,i ), mortality rates of reserve biomass (m R,i ) and one of the shape parameters of the conversion rate from reserve to green biomass (c i ).Black points denote the non robust parameter sets with 0.1≤P 3,ι ≤0.5, and red filled circles denote the robust parameter sets with P 3,ι >0.5.The axes show the entire parameter space that was sampled in model C. The clustering of z r,i , w R,i , and c i is caused by the constraints in parameter space and the plausibility check (see Sect. 2.6).coexistence again.In this case, coexistence results from the combination of niche differentiation and environmental disturbance.The latter fits in the context of the Intermediate Disturbance Hypothesis (Connell, 1978;D'Odorico et al., 2008;Grime, 1973;Huston, 1979).The species specific flood resistance in model C allows for ecological differences in the response to disturbance and outbalances too strong advantages from the differences in the phenology, and thus enhances coexistence (Roxburgh et al. 2004).
Although both models differ in their structure and coexistence mechanisms, the ensemble statistics of mean hydrologic variables like transpiration and depth to ground water are surprisingly similar between models A and C (Fig. 6).This is owed to the fact that the hydrological model is the same in both A and C.However, the differences between the two models become apparent, when considering the variation in the time series for both hydrological and ecological variables (depth to ground water, green and reserve biomass) of the system and its sensitivity to environmental change (here: change of the Hurst exponent).The more complex model C shows higher variation in the variables, and is more sensitive to environmental change than model A. This is a logical consequence of the modelled co-existence mechanism.In model C, the flood has both indirect (via the hydrosystem as resource) and direct (as disturbance) impacts on the plant species.Thus, both reserve and green biomass of the different species are independently linked to the flood fluctuations.As a result, species abundances change over time, sometimes with a prevalence of the water conserving species, sometimes with prevalence of the water demanding species.Thus, transpiration and the resulting ground water level vary accordingly.In model A, however, the flood influences the ecosystem merely via the reserve biomass (no direct impacts on the green biomass).The reserve biomass is able to act as a buffer and to stabilize the entire system (green biomass, ground water depth).
The results on the influence of the Hurst exponent also give rise to some conclusions on the adequateness of the Intermediate Disturbance Hypothesis (IDH) in ecohydrological systems along ephemeral rivers.The IDH primarily argues with the frequency of the disturbance.Our results indicate, however, that the autocorrelation in the varying water supply and so the duration of related disturbance events (cumulative water stress during dry periods, repeated disruptive floods) are crucial for the impact on species coexistence and resilience as well.In this ephemeral ecosystem, considering solely the frequency would reach too short.The importance of autocorrelation/red noise has also been shown in the context of species survival.Schwager et al. (2006) for instance, showed that autocorrelation can be stabilizing or destabilizing depending on the species' ecological traits.
The results of our study suggest that the assumptions on the functional traits of the species in the plant communities (e.g.regarding resource utilization, flood resistance) and so on the mechanisms of competition/coexistence can influence the modelled hydrology.Furthermore, we find hints that the distribution of mean hydrologic variables in this system is probably driven by the applied hydrological model, whereas the distribution of fluctuations (here: coefficient of variation) is probably driven by the assumed ecological interactions.
Our forward simulations with different Hurst exponents show that not only the stochasticity of the environmental disturbance (the flood) influences the coexistence of the three species, but also the cyclicity of periods with high and low floods (long term memory, see Sect.2.3) plays an important role.Most hydrological processes are characterized by long term memory processes (Montanari et al., 1997), which lead particularly in arid regions to extended periods of unusually small or strong events ("Joseph Effect", Mandelbrot and Wallis, 1968).Our model suggests that this hydrologic characteristic might have important influence on ecosystem structure.This finding is in line with other results showing that the fine structure of environmental fluctuations can alter systems dynamics, qualitative trends or ranking orders among scenarios with serious implications for management (Frank, 2005;Schwager et al., 2006).Furthermore, the two models A and C show differences in the sensitivity of species coexistence against a change in the Hurst exponent.While model C reveals a strong sensitivity and a loss of coexistence, model A is found to be rather robust.The reason for this difference is again the buffer capacity of the reserve biomass in absence (model A) and presence (model C) of direct disturbance effects of the flood on the plant species.The two models A and C can also be interpreted as two types of plant communities which differ in the impact of floods on their species (e.g.indirect only; indirect and direct).But note that both models that successfully modelled robust coexistence are still abstract representations of ecological and hydrological processes along ephemeral rivers.Thus, only limited knowledge of the actual mechanisms is implemented.Such generic models that focus on essential aspects are known to be crucial for integration and analysing consequences of feedback loops when entering new interdisciplinary fields (Baumgartner et al., 2008).This allows formulating new hypotheses, which can then be tested by more complex and structurally realistic models.In our context, additional intra-or interspecific effects (like age dependent rooting depth) might be active in maintaining the observed coexistence.Potentially, a lot more mechanisms can enhance the three species coexistence like random individual effects or multi dimensional tradeoffs (Clark et al., 2007).Thus, our models are just a subset of possible abstraction, which might all reproduce the observed coexistence.In fact it might be impossible to find the "right" model.The coexistence constraint did not limit the possible parameter space enough to lead to a unique ecohydrologic response.However, our models shed light on possible options.They also give hints towards which variables could be measured to increase the understanding about the involved mechanisms.

Conclusions
The modelling of three species coexistence in a water limited environment is challenging because feedbacks between ecology and hydrology have to be implemented in an appropriate way.The present study introduced a model that facilitates the investigation of effects of model structure and parameter uncertainty on ecology and hydrology of the water limited system along ephemeral rivers.We applied a range of model versions with a varying degree of included information.Given that only two of three models led to robust three species coexistence, we conclude that the driving coexistence mechanism is defined by the model structure.On the other hand, the robustness check of the parameter sets leading to three species coexistence indicates that the success of the underlying coexistence mechanism is controlled by the combination of the population parameters.Further, depending on the model structure the flood can act as water resource or environmental disturbance or a combination of both.When acting as environmental disturbance the change in long term memory strongly affected the robustness of the parameter sets.Therefore, we conclude that the long term memory of hydrological processes is important in water limited ecosystems.In this study, we applied the same hydrological concept for all model versions and only changed the complexity of the ecological model.Considering that the distribution of average values of transpiration and ground water table were similar but not their distribution of fluctuations, we conclude that the ensemble statistics of average values of hydrologic variables are probably influenced by the applied hydrological model, whereas the ensemble statistics of fluctuations of both are probably controlled by the applied ecological model.
Our study shows that the species composition in the plant community strongly influences the stability properties of the ecohydrological system (e.g.variation in transpiration and ground water depth; variation in reserve and green biomass; sensitivity of species coexistence to change in the Hurst exponent).This stresses the necessity to consider explicitly species composition and functional interactions in the ecosystem when assessing the impact of climate or land use change on water resources and vegetation along ephemeral rivers.This is particularly important in systems where the floods have direct destructive impacts on the vegetation.Here, models are essential that explicitly take into account such disturbance effects (such as model C).The relative importance of the species composition for understanding eco-

Fig. 2 .
Fig. 2. Water balance of an ephemeral river with shallow aquifer.The intermediate zone denotes the layer where saturated and unsaturated conditions alternate frequently.The arrows denote the transpirational demand for each species T WS,i (Eq.11), the infiltration flux Q I , and the ground water recharge GWR (Eq.6).

Fig. 3 .
Fig. 3. (a) Flood volume and (b) duration at gauging station Schlesien from 1981 to 2006.Data are provided by the Department of Water Affairs (DWA) in Windhoek.

Fig. 5 .
Fig. 5. Results of the parameter sampling.Histograms show the relative frequency of parameter sets resulting in P 3 ≥0.1 with model C and A (H =0.75 for both).

Fig. 6 .
Fig. 6.Relative frequencies of ensemble mean total transpiration (left column) (Eq.31) and ensemble mean depth to ground water (right column) of parameter sets with P 3 ≥0.1 for model C (upper row) and model A (lower row).

Fig. 7 .
Fig. 7. Relative frequencies of ensemble CV of total transpiration (left column) (Eq.33) and ensemble CV of depth to ground water (right column) of parameter sets with P 3 ≥0.1 for model C (upper row) and model A (lower row).

Fig. 8 .Fig. 9 .
Fig. 8. Medians of Ttotal (left column) and zunsat (right column) of parameter sets with 0.1≤P 3,ι ≤1.The bin size of x-axis is 0.01.Results of model C are shown in the upper row and model A in the lower row.The linear correlation coefficients are (line by line): 0.34, −0.38, −0.08, −0.06

Fig. 11 .
Fig. 11.Typical time series of model A: (a) reserve and green biomass, (b) flood volume and depth to ground water.

Fig. 12 .
Fig. 12.Typical time series of model C: (a) reserve and green biomass, (b) flood volume and depth to ground water.

Table 1 .
Ecology of the three main tree species along the middle part of the Kuiseb River.

Table 2 .
Hydrological parameters (for soil and flood shape) used for all model versions (Table

Table 3 .
Transpiration rates for each species.

Table 4 .
Levels of complexity in model A-C.Phenology and flood resistance can be implemented species specific or same for all species.

Table 5 .
Ecological parameters that were calibrated and their range.

Table 7 .
Symbols used in this study, i denotes the reference to a species.

Table 7 .
Continued.Water volume in the alluvium where no water is available for roots of species i m 3 ha −1 27 V Stress,i Water volume in the alluvium that leads to water stress of species i