Parameterization and uncertainty in coupled ecohydrological models

In this paper we develop and apply a conceptual ecohydrological model to investigate the e ﬀ ects of model structure and parameter uncertainty on the prediction 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 5 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 ﬂood. The ecosystem is modelled with a population model, and repre-sents three dominating riparian plant populations. In appreciation of uncertainty about population dynamics, we applied three model versions with increasing complexity. Pop- 10 ulation 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 di ﬀ erent coexistence mechanisms, and reacted di ﬀ erently to change of long term

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. Pop-10 ulation 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

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 understand-10 ing 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 15 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 20 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 25 Gross, 1992; Grimm and Railsback, 2005), representing a bottom-up 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 4157 to be analysed. The respective models, however, are often complex that makes parameterization 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. There-5 fore, we follow a top-down 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 at the individual scale (Moorcroft, 2003;Frank and Wissel, 2002;Fahse et al., 1998;Heinz et al., 2005;Hanski, 10 2002, 2004). Most ecohydrological models work at the population scale Porporato et al., 2001;Camporeale and Ridolfi, 2006;Ridolfi et al., 2000). Direct parameterization 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 15 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 parameterized.
In this study, we address this parameterization problem by using pattern-oriented 20 modelling, in that we adjust species parameters such that the resulting model reproduces the observed coexistence. Models have been parameterized based on information of presence or absence of plant species before. When the existence criterion will be extended to several species it is called coexistence, and also observed coexistence has been used to evaluate (at least qualitatively) the validity of ecohydrological 25 models (Laio et al., 2001;Rodriguez-Iturbe et al., 1999). 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). It is required, that mechanisms supporting coexistence are related through precise order of parameter combinations. A given model only allows for coexistence, if its parameters meet certain conditions. 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 5 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. However, many hydrologic time series are characterized by auto-correlated and longterm-memory processes (Montanari et al., 1997;Hurst, 1951), particularly in arid environments. This 10 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 15 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 20 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 inter-25 mittent floods. At the same time, strong floods lead to uprooting of riparian vegetation and increased mortality. There is no 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 4159 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 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, 5 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 mean hydrologic variables (groundwater level and transpiration) were similar in all models, their variability depended both on the model structure and the 10 parameters sets. This points at the difficulty to parameterize an ecohydrological model in real world applications. However, our model gives clear indications, what measurements are most effective for improving the necessary process understanding. 15 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 20 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 25 creates a living environment for riparian vegetation. The flood is influenced by up-4160 stream farm damns and the ground water table is influenced both by plants and human consumption.

Ecosystem
Vegetation around the river channel consists to 80% of only three coexisting species: Camel Thorn (Acacia erioloba), Ana Tree (Faidherbia albida) and Wild Tamarix 5 (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 10 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) (Fig. 2). Its thickness and width vary along the river. The alluvial aquifer is 20 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 infiltra- 25 tion. The dynamics of flood water infiltration were investigated by Dahan et al. (2008).

4161
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 fluxes are very similar although each flood is characterized by different flow conditions. 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 .

10
The water balance for this segment is written as Eq. (1) where ∆W S(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 15 S unsat (t) = A seg · z unsat (t) · θ(t) for the unsaturated storage, (2a) 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.
The change in unsaturated storage was calculated as with I(t) denoting the infiltration, GW R(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 5 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)=2400 m 3 /d ha. Therefore, the infiltration depends only on flood duration D(t) (Eq. 20) and the specific infiltration flux Q I (t): 10 Flood durations is calculated as a function of flood volume in Eq. (20) (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 15 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): For plants where the roots reach the groundwater, transpiration originates from both 20 the unsaturated and the saturated zone. The unsaturated part is calculated as

4163
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), T WS,i (t) the transpirational demand for each species Eq. (10), 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 by plant roots. The water in the unsaturated storage available for transpiration of species i depends on its rooting depth z r,i (t): 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): (10) 15 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 ∆S 20 where Q GW (t) is the ground water flow and T GW (t) the transpiration of all species from ground water Eq. (14). 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: 5 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 15 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 aver-20 age (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 4165 that are present in many hydrologic processes (Montanari et al., 1997;Hurst, 1951). 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 : 10 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: 15 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 20 series were generated with FARIMA(p=1, d =0. 25, q=1) and adjusted to the observed mean annual flood volume µ Flood = 3 269 000m 3 , and thus yielding V Flood (t) = e (FARIMA(1,0.25,1)+log(µ Flood )) .
4166 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, 15 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 20 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. (24), and ε i (t) is the unitless water stress function Eq. (25), 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 5 available amount of water (Eqs. 24 and 25). The first term of Eq. (21), (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 from the existing reserve biomass.
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 phe-10 nology. In the first case, Eq. (21) 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 l s i is the unitless leaf shedding factor and ranges from 0 to 1 of species i . l s i =0 15 corresponds to no leaf shedding at all, and 1 to complete leaf shed. Usually, leaf shed is not complete, so l s i takes a value between 0 and 1. The formation of reserve biomass takes place at the end of each season t: where fr i (t) is the unitless flood resistance of species i and ranges from 0 to 1 (Eq. 26), 20 see below). It denotes the vulnerability of a given species to being uprooted and washed away by a flood of given magnitude. fr i (t)=0 corresponds to complete removal of reserve biomass by the flood. In the dry season, fr i (t) is set to 1. The parameter m R,i denotes the mortality of the reserve biomass, and w R,i the rate of conversion from green into reserve biomass. Both are constant over time and unitless. The first , 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 5 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. 23).
The three parameters conversion rate w G,i (t), water stress ε i (t), and flood resistance 10 fr 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 WS,i (t)) (see Sect. 2.2, Eq. 8) 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.

20
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 species-specific because it depends on the species root depth Eq. (9). V Stress,i is also a species-specific parameter: the lower V Stress,i the more drought tolerant is this species.

5
The flood resistance fr 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. (23). We model it as a linear function of the flood volume (V Flood with unit m 3 /ha), which was generated by Eq. (19).
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 15 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 20 representation of the ecosystem. Since little is known about the ecological parameters in the Kuiseb River, any model would be comparatively simple. We compare 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 Tree (Table 1). For this we added two parameters (l s Cam and l s Ana ), which increased the degree of complexity (model type B). Finally, in model C, we included more knowledge regarding the difference in flood resistance between species, thus allowing the parameter fr i to be species specific.
In summary, model type C included the most ecological information, strongest con-5 straints 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 parameterized the models accordingly, as pointed out in the next section.

Parameter sampling
10 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 parame- 15 ters 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 den-20 sity 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 The sampling procedure is illustrated in Fig. 4. For each sampled parameter set (Ω ι ) 25 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 4171 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 parameter set Ω ι only led to coexistence under very specific flood conditions, 5 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 10 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 15 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 20 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 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 15 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 4173 term memory of the flood generation algorithm, by decreasing and increasing the Hurst exponent (Eq. 18). 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). Table 6 shows in per cent how many of the 150 000 sampled parameter sets led to 5 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).

10
Model B was not subject to further investigations because of the low number of parameter sets leading to three species coexistence. 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 15 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 ensem-20 ble 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 25 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 per-centile 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 ensem-5 ble averages of transpiration and depth to ground water, they were much different with regard to 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).

10
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 15 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 exists 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 exists between P 3,ι and the CV of 20 transpiration (r 2 =−0.27) and ground water table (r 2 =0.14). No such relation exists 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 25 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 4175 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 5 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 relative 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. 15 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 infor-20 mation, 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, competition models from ecology (e.g. Lotka, 1925;Volterra, 1926) underpinning that different species can only coexist under the precondition that inter-specific competition (between species) is weaker than intra-specific competition (within the same species). As a result, coexistence can only be found for appropriate parameter combinations that 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 coexistence. In this study, only models A and C allow for 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 15 influencing green biomass and, hence, transpiration demand of the individual species (see Eq. 10). 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) phenol-20 ogy. 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 25 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. In model C, another coexistence mechanism is enabled, only by allowing for species 4177 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, 5 hence, can mediate 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 (D'Odorico et al., 2008;Huston, 1979;Grime, 1973;Connell, 1978). The species specific flood resistance in model C allows for ecological differences in the response to disturbance and outbalances too 10 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 15 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 20 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 25 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 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) 15 and so on the mechanisms of competition/coexistence can influence the modelled hydrology. Furthermore, ensemble statistics of mean hydrologic variables in this system are driven by the applied hydrological model, whereas the ensemble statistics of fluctuations (here: coefficient of variation) is driven by the assumed ecological interactions.
Our forward simulations with different Hurst exponents show that not only the 20 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

25
Effect", Mandelbrot and Wallis, 1968 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 5 (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 coexistence are still abstract representations of ecological and hydrological processes along ephemeral rivers.

10
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 15 intra-or interspecific effects 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 20 co-existence 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.

25
The modelling of three species coexistence in a water limited environment is a 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 the observed three species coexistence, we conclude that the 5 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 10 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 average 15 values of transpiration and ground water table were similar but not their fluctuations, we conclude that the average values of hydrologic variables are mainly controlled by the applied hydrological model, whereas the fluctuations of both are controlled by the applied ecological model.
Our study shows that the species composition in the plant community strongly in-20 fluences 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 25 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 ecohydrological sys-4181 tems, however, came only to light through the subsequent process of changing the model structure and comparing their outcomes.

5
The number of adult individuals in population i (N Ind,i ) was calculated to define the coexistence criterion for the parameter sampling (Sect. 2.6): where R 1,i denotes the reserve biomass of one adult individual of population i . We simplified the shape of an individual (reserve biomass above and below subsurface) to 10 be a right circular cylinder with maximal trunk radius r i , maximal height h max and wood density ρ i :  (Curtis and Mannheimer, 2005), b (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) 4187 Table 2. Hydrological soil and FARIMA parameters used in this study. 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.    Ground water inflow 14.9 m 3 ha −1 season −1 (dry) 12 20.9 m 3 ha −1 season −1 (rainy) Q GW Ground water flow m 3 ha −1 season −1 12 Q L Lateral ground water inflow 869.2 m 3 ha −1 season −1 12 4193         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).