A model study of terraced riverbeds as novel ecosystems

Riverbed terracing has been introduced in ancient times to retain water and soil, to reduce hydrological connectivity and erosion and to increase primary and secondary productivity of agro-ecological systems. These presently abandoned human-made landscapes have become novel ecosystems and a potential source of ecosystem services to humans in drylands. We use a mathematical-modeling approach to study factors that regulate terraced riverbeds and affect community and ecosystem attributes such as productivity, functional diversity and resilience to droughts. We introduce a model that captures the relationships between rainfall pattern, runoff coupling between adjacent terraces, and vegetation growth, taking into account competition for water and light. We found that a large number of weak rainfall events results in lower total biomass and functional diversity across the terraced riverbed compared with a few strong rainfall events. We further analyzed the filtering of species traits from pools of functional groups that make different tradeoffs between investment in above-ground biomass to capture canopy resources and investment in below-ground biomass to capture soil resources. Pools characterized by concave tradeoffs give rise to higher functional diversity, lower biomass production and lower resilience to droughts, as compared with convex pools. New empirical studies are needed to test these model predictions.

gradually disappeared when farming practices were abandoned, and the natural buildup of drainage networks has increased soil erosion and transport of sediment by runoff. These processes are especially visible in semi-arid environments with sparse plant cover that were subjected to irregular rainfall events and long periods of drought. The net effect of the soil erosion processes is a significant reduction in soil quality and its function in C and N cycling.
The negative effects of terrace abandonment are partially compensated by natural ecological processes that form the present ecological legacy of the terraced system. Abandonment is followed by re-establishment of the dry valley through natural vegetation that grows on the watershed slopes. Annual and perennial plants species and their associated animal communities create new ecosystems in conjunction with the abiotic conditions of the terraces dictated by the hydrological and geomorphological legacies, termed novel ecosystems 6 . The self-organized ecosystem offers an increase of water retention and decrease in soil erosion mainly by woody plants that function as hydrological engineers and enhance the soil infiltration capacity. In addition, herbaceous and woody plant species improved soil structure and function due to increased organic content resulting from litter fall and decomposition. Thus, abandonment of farming practices, followed by the re-colonization of natural vegetation, has improved soil functions in relation to water and nutrients fluxes, and created effective protection against soil erosion 12,13 .
Terraces, along with their hydro-geo-eco legacies, are crucial components of the ecological landscape of drylands, but very little is known about the soil-vegetation feedbacks underlying hydro-geo-eco processes of terraced areas and their effect on productivity and diversity. Despite their long periods of abandonment, ancient terraces constitute a reusable landscape capital for combating desertification and mitigation of climate change by enabling local retention of water and organic matter and sustaining plant communities 8 . Restoring their function, re-designing their water-connectivity architecture, and re-managing them in accord with rainfall variability, can result in productive and resilient dryland ecosystems.
Assessing the potential of ancient terraces as a landscape capital can highly benefit from mathematical modeling and model studies. Models of dryland ecosystems [14][15][16] have already proven to be instrumental in understanding self-organized vegetation patchiness [17][18][19] and functional diversity of plant communities 20 . The aim of this paper is twofold. We first use a meta-ecosystem modeling approach [21][22][23] , and earlier dryland-ecosystem models 24 , to introduce a new terraced-riverbed model that captures the relationships between rainfall pattern, runoff coupling between adjacent terraces, and the differential plant-community dynamics in the terraces. We then use the model to study whole riverbed responses, in terms of biomass production and functional diversity, to different rainfall regimes, varying the tradeoff form that characterizes the species pool. The specific model we study is motivated by the structure of the extensive terraced system in the Negev desert ( Fig. 1), wherein an area of about

A Meta-ecosystem Model for terraced Riverbeds
Our starting point is a model for a single terrace (ecosystem) that describes a water-limited plant community 24,26 . It consists of two water variables, below-ground or soil water, W, and above-ground or surface water, H, and biomass variables B B B , , , m 1 2 … , representing the above-ground biomass of m functional groups. We focus on functional groups of species that make different compromises in their relative investments in above-ground biomass (shoot), to capture canopy resources, and below-ground biomass (root), to capture soil resources. Specifically, we focus on two functional traits, the maximal shoot biomass k, and the root-to-shoot ratio E, and introduce a dimensionless tradeoff parameter χ ∈ [0, 1], defined implicitly through the relations: min m ax min where K min , K max and E min , E max determine the ranges of the two traits in the pool of functional groups, and  is a positive constant (Nathan et al., 2016). Thus, 0 χ = represents a functional group that invests mostly in shoots to capture light, , while 1 χ = represents a functional group that invests mostly in roots to capture soil water, The distributions of all other functional groups χ i , in between these two extremal groups, strongly depend on the value of  as the tradeoff curves in Fig. 2 show.
The significance of  can be understood by considering an ideal functional group of species K E ( , ) max m ax , a "Darwinian demon" 27 , that benefits from high values of both K and E (red point in Fig. 2). Functional groups that are closer to this ideal group are expected to have competitive advantage over other groups. The closest functional group can be calculated by minimizing the distance function max m ax 2 2 with respect to χ, where χ K( ) and E( ) χ are given by a (1 ) and b (1 ). For  = 1 (linear tradeoff), such a calculation gives the intermediate value www.nature.com/scientificreports www.nature.com/scientificreports/ for the closest, and thus most competitive, functional group. As  is decreased below unity (convex tradeoff) a distinguished set of functional groups around mid χ becomes significantly closer to the ideal group K E ( , ) max m ax than all other groups. By contrast, as  is increased beyond unity (concave tradeoff) the distance to the ideal functional group becomes more uniform across the pool, implying a wider range of functional groups with nearly equal competitive capabilities. Note that for sufficiently high  values the shortest distance to the ideal group occurs at the two extreme groups, χ = 0 and 1 χ = .
The biomass variable associated with the ith functional group is B B( ) , that is the biomass obtained by solving the model equations with the particular parameter values K K( ) , keeping all other parameters fixed. For simplicity, we assume that in each terrace all water and biomass variables are spatially uniform. This assumption holds for homogeneous systems with weak pattern-forming feedbacks 19,28 .
The terraced riverbed (meta-ecosystem) to be studied here is obtained by coupling the single-terrace models for adjacent terraces through the above-ground water variable to account for runoff contributions from upslope terraces and loses to downslope terraces. Note that we did not include sediment transport from upslope terraces as this is a secondary effect. The meta-ecosystem model for a riverbed with n terraces reads where the index i runs over the m functional groups, the index j runs over the n terraces, B ij is the above-ground biomass of the i th functional group in the j th terrace, and W j and H j represent, respectively, the contents of soil water and of surface water in the j th terrace. The factor ij Λ in equation a (1 ) is given by: and represents growth attenuation of species that suffer from competition for light (Nathan et al., 2016). This competition is strongly affected by the maximal-shoot biomass parameters, K i , which affect the late-stage growth of the plant species. These parameters represent species-specific growth limitation, assumed here as resulting from self-shading, and thus related to specific leaf area. Competition for water is captured by the water-dependent biomass growth factors ( ) and by the water-uptake rates Γ ∑ + ( ) . That competition is strongly affected by the root-to-shoot ratios, E i , of the different functional groups. The source of soil water W j in any terrace comes from the infiltration of surface water H j at a rate I. The source of surface water comes from precipitation at a time-dependent rate P and runoff contributions from adjacent upper terraces and from surrounding hillsides. The latter contribution often comes predominantly from the hillslopes that surround the uppermost terrace = j n and is modeled in Eq. (2d) by the term P α , where α is a function of the precipitation rate given by This form captures the diminishingly small runoff contribution at low P when the soil is dry, α α ≈ P P P 1 3  , the surge in that contribution as P increases, and its attenuation at high P, as the soil becomes saturated and the runoff-contribution approaches a linear dependence, α ≈ in Eq. (2c), where the parameter D provides a measure for the coupling strength. An illustration of the factors that affect surface-water dynamics in a terraced riverbed is shown in Fig. 3.
The reader is referred to Table 1 for a list of all model parameters, their meanings, units and numerical values. A few more comments are in order. We assumed that the biomass loss parameter M i , which represents mortality and additional factors such as grazing, is uniform across the community M M i = . This is unlike an earlier study 20 where a distinction between short and long grazing histories was made. We have not modeled transport of soil water between the terraces as it is much slower that the transport of surface water 8 . We have not included an evaporation term in the surface-water Eq. (2c,d) because of the fast conversion of surface water to soil water through infiltration. For simplicity, the runoff rate D is considered to be a constant. In practice, runoff is slowed www.nature.com/scientificreports www.nature.com/scientificreports/ down by vegetation and D is biomass dependent. We note that D is related to the uniform extension = − + l x x j j 1 of the terraces and to their uniform height d z z

Results
Biomass distributions along the tradeoff axis. Integrating the model equations for sufficiently long times we obtained asymptotic biomass distributions, as Fig. 4a shows, for a riverbed consisting of N 9 = terraces with a linear tradeoff ( = 1  ). The asymptotic pulse-shape distributions contain information about three community-level properties: the total biomass (pulse area), functional diversity (pulse width) and community composition (pulse position) 20 . Inspecting the distributions, from the uppermost terrace j 9 = to terraces downstream, we reveal declines in total biomass and functional diversity, which occur mostly across the first three terraces, and a composition shift to functional groups that invest more in below-ground biomass (higher χ). Figure 3. An illustration of surface-water dynamics in terraced riverbed. Gain of surface water in the jth terrace is due to precipitation, P, and runoff contribution, + DH j 1 , from the adjacent upper terrace. Loss of surface water in the jth terrace is caused by infiltration into the soil, IH j , and by runoff, DH j , to the adjacent lower terrace. The uppermost terrace receives a runoff contribution, αP, from the surrounding hillslopes. , as Fig. 4b shows for the three uppermost terraces (results for the lower terraces are not shown as the variations become negligibly small). The nonlinearity, however, introduces two additional effects. As the tradeoff changes from convex to concave the change in community composition (pulse position) along the riverbed increases and the total biomass decreases. This is in line with the view of a concave tradeoff as representing a pool with a wide range of equally competitive functional groups that are farther away (in trait space) from the ideal functional group (red circle in Fig. 2), as compared to a convex pool.
The composition change of the community along the terraces implies higher functional diversity of the whole terraced riverbed as compared with any individual terrace. This effect is particularly significant for concave pools for which the composition change downstream is larger.
Thus, riverbeds with concave pools of functional groups are expected to show a significantly higher functional diversity, in comparison to convex pools, but lower total biomass.
Response to droughts. The lower biomass associated with concave pools may bear on the response of the community to droughts. Figure 4c shows the responses to 10-years long and 60-years long periods of dry climate (MAP reduction from 100 mm/y to 40 mm/y) applied to concave and convex pools. While convex pools show a moderate biomass decline for both periods of dry climate, concave pools show a severe decline for the longer dry period. Continuation of the model simulations for a 10 years period of normal climate (100 mm/y) after the 60 years of dry climate show very fast recovery for convex pool (100% of the original biomass) and very slow recovery (8.4% of the original biomass) for concave pool. ( 1). The width of the pulse (functional diversity) and the area it encloses (total biomass) decrease as we go downstream from the top terrace j ( 9) = , and the pulse position shifts to higher χ, reflecting a compositional change towards higher investment in below-ground biomass at the expense of above-ground biomass. These changes become negligible below the = j 6 terrace, indicating weak water-flow connectivity between the lower terraces.  Table 1  www.nature.com/scientificreports www.nature.com/scientificreports/ Changing the number of rainfall events keeping MAP constant may have strong effects on functional diversity and total biomass, as Fig. 5 shows. The functional diversity and total biomass calculated for the entire riverbed sharply decline as the number of rain events increases and their intensity decreases, and approach constant values as the rainfall drops to levels for which runoff contributions from the surrounding hillslopes of the terrace become negligible. These behaviors hold for all pool types, convex, linear and concave. Comparing the functional diversity and total biomass for these pool types, at any given number of rain events, we recover the opposite trends that functional diversity and total biomass follow when the tradeoff-nonlinearity  is changed, as discussed earlier.

Discussion
Studies of terraced riverbeds have focused mainly on their function as agricultural systems and on the rehabilitation and conservation of abandoned riverbeds to prevent soil erosion and fertility loss 29 . Only a few field studies have focused on the ecological aspects of plant colonization in abandoned terraced riverbeds and the novel ecosystems this colonization has created 30 . In this paper we used a mathematical-modeling approach to study factors that regulate these novel ecosystems and affect ecosystem attributes such as productivity and diversity. We applied a meta-ecosystem modeling approach 21,22 to earlier dryland-ecosystem models 24 , to obtain a model for a terraced riverbed that captures the relationships between rainfall pattern, runoff coupling between adjacent terraces, and the differential plant-community responses in the terraces in terms of biomass production and functional diversity. Our model studies show that rainfall patterns and functional properties of species pools regulate biomass production, functional diversity, and responses to extreme droughts, at the levels of a single terrace and whole riverbed, as discussed below 31 . Based on our model studies, rainfall and its redistribution by runoff are the main drivers of ecosystem development in terraced dry riverbeds. The rates of the cascading process of water flow, rainfall → runoff → soil moisture → vegetation, are regulated by the rainfall pattern, i.e. by the temporal distribution of rainfall events in terms of frequency and magnitude. Varying the number of rainfall events in a year, keeping the mean annual rainfall (MAP) constant, we found that a large number of weak rainfall events results in low total biomass and low functional diversity (Fig. 5). We attribute these outcomes to the negligible hill-slope runoff that is generated by weak rainfall events. The low total biomass is a consequence of the low soil-water content in the absence of runoff contributions. The low functional diversity is a consequence of the high similarity among terraces; as the main water resource in all terraces is direct rainfall, their community structure is similar too. By contrast, in the case of a few strong rainfall events, the generated runoff provides an additional source of water that increases not only the total biomass but also the functional diversity. This is because of the decreasing runoff contributions to terraces down the riverbed, and the consequent different communities that establish there: functional groups specializing in capturing light (low χ) in the water-rich upper terraces vs. functional groups specializing in capturing soil-water (high χ) in the water-poor lower terraces (Fig. 4). The functional diversity at the riverbed level takes into account the low χ functional groups in the uppermost terrace and the high χ functional groups in the lowermost terrace. It is, therefore, necessarily higher than the functional diversity in the case of many weak rainfall events, where similar communities establish in all terraces. We thus conclude that fewer stronger rainfall events result in more functional ecosystems. , and concave > ( 1)  tradeoffs. Functional diversity (a) and total biomass (b) both decrease as the number of rain events increase (keeping MAP constant). For a given number of rain events functional diversity decreases (a) while total biomass increases (b) as the tradeoff changes from concave (blue) to convex (red).
We are not aware of any experimental study aiming at testing directly the relationship between rainfall pattern and properties of properties of terraced-riverbed ecosystems. However, some information about this relationship can be inferred from studies of human-made systems termed "Limans" in the central Negev, which can be regarded as systems consisting of a single terrace per watershed [32][33][34] . In particular, it has been found 35 that most annual rainfall events are in low magnitude and high frequency, and that the low-magnitude high-frequency events redistribute the rain within the slopes of the drainage area only. Runoff contributions to the riverbed occur only during rare high-magnitude events. Thus, rainfall events regulate ecological responses, such as species diversity, mostly on slopes, while the regulation within the riverbed occurs only during rare intense rainfall events. These results are in line with the model predictions shown in Fig. 5, namely, that biomass production and functional diversity are hardly affected by the rainfall pattern in the high-frequency range, but show substantial changes in the low-frequency range. We adopted in this study a trait-based modeling approach, which directly bears on functional diversity and ecosystem function [36][37][38] . The approach allowed us to study the filtering of species traits from a set of species pools and the local community structures that emerge from interspecific competition for water and light under various rainfall regimes. The species pools we consider consist of functional groups that make different tradeoffs between two resource acquisition strategies: investment in above-ground biomass to capture canopy resources, and investment in below-ground biomass to capture soil resources. Two major forms of such tradeoffs, and thus of species pools, can be distinguished according to their shapes in trait plane: convex tradeoffs ( 1  < in Fig. 2) and concave tradeoffs ( 1  > in Fig. 2). Empirical information about tradeoff shapes in plant ecology is very limited 39 , but a recent study provides insightful information for the particular tradeoff we consider here 40 . This study focused on specific leaf area (SLA) and specific root length (SRL) as functional traits associated with above-ground and below-ground resource capture, and while graphs of SLA vs. SRL suggest both positive and negative relations between the two functional traits, replotting the data as graphs of SLA vs. SRL/SLA, to better correspond to the model's functional traits K and E, results in a clear concave tradeoff, (see Fig. A2 SM). This result supports the tradeoff relation (1) that we assume between K and E, and points towards the possible existence of concave tradeoffs in dryland species pools.
As pointed out earlier, a heuristic view of the difference between convex and concave pools is provided in terms of the distance χ s( ) in the E K , trait plane between a functional group represented by a point χ on the tradeoff curve and the "ideal functional group" E K ( , ) max max , which does not suffer from making a tradeoff [see Eq. (1c) and red dot in Fig. 2]. Short distances imply superior environmental-filtering and competitive capabilities, which suggests that convex pools provide advantages to small sets of intermediate functional groups whose distances to the ideal functional group are significantly shorter than those of the extremal groups (i.e. near 0 χ = and 1 χ = ). By contrast, moderately concave pools give rise to higher similarity among functional groups, as the distance s( ) χ of intermediate groups is now comparable to that of the extremal groups. As a consequence, we may expect functional diversity to increase and environmental filtering to decrease as compared to convex pools.
These expectations and their implied consequences for ecosystem function are supported by our model results. As Fig. 5 shows, concave pools give rise to higher functional diversity (Fig. 5a) and lower biomass production (Fig. 5b) as compared with convex pools, for all rainfall patterns considered. These trends appear already at the single terrace level (Fig. 4b), but are more pronounced at the whole riverbed level because of the decreasing runoff contributions to lower terraces and the consequent shift in community structure to functional groups of higher χ. The shift itself increases the whole-riverbed functional diversity, and the specific shift to functional groups that invest less in biomass production contributes to lower productivity. The better performance of convex pools in terms of biomass production is also reflected in the response to periods of dry climate. As Fig. 4c shows, a long period (60 y) of reduced MAP has a moderate effect for a convex pool but dramatically reduces the biomass in the case of a concave pool. In addition, the recovery time of a concave pool, when the original MAP is resumed, is much longer than that of a convex pool.
We have studied here moderately nonlinear tradeoff curves, focusing on the difference between convex and concave pools. For such tradeoffs we found intermediate community structures that exclude functional groups with extreme investment either in shoot or in root (Fig. 4). However, tradeoffs in practice may assume stronger nonlinear forms. Consider, in particular, a highly concave tradeoff, as shown in Fig. 1 with  =4 or as the empirical data in Fig. A2 of the SM shows. The simple heuristic reasoning based on the distance function χ s( ) suggests, in this case, the possible split of the community into two extreme sub-communities. The first sub-community consists of functional groups with extremely low investment in root ( ≈ E E min ) and variable investment in shoot, whereas the second sub-community consists of functional groups with extremely low investment in shoot (K K min ≈ ) and variable investment in root. The former may be expected to dominate at relatively high MAP, whereas the latter at a relatively low MAP. At intermediate MAP, bistability of the two sub-communities may possibly exist, implying abrupt shifts in community structure as MAP decreases or increases.
Prevailing trait-based frameworks 38 focus primarily on community assembly from a fixed species pool, paying little attention to species pools that differ in the nonlinear forms of the tradeoffs that characterize them. Our results show that convex and concave tradeoffs result in significantly different community-level properties. These essential differences call for additional theoretical and empirical studies of convex vs. concave pools, in terraced riverbeds and in other dryland meta-ecosystems, and their implications to functional diversity, biomass production, resilience to droughts, and ecosystem management. To this end, the modeling approach pursued here in the context of terraced riverbeds can be modified or extended in several directions, including the consideration of additional tradeoffs, such as fast growth versus ability to withstand a resource stress 41 , and the modification of the water-connectivity architecture.
Hydro-geo-eco legacies shape the long-term structure, composition, and function of terraced systems. Using a model that captures essential aspects of these legacies we uncovered basic determinants of plant-community structure, such as the relation between number and intensity of rainfall events, and the nonlinear forms of tradeoff relations. Mechanistic insights obtainable in this way, can increase our scientific-based knowledge for effective management of terraced riverbeds. We suggest that various ecosystem services, such as regulation of water flow and primary production, can be achieved by re-maintenance of these cultural landscapes according to these insights.

Methods
We study the model Eq. (2) numerically using precipitation regimes that mimic a few essential aspects of typical intra-annual rainfall patterns in the Sde-Boker area of the Negev Desert, Israel. That pattern, shown in Fig. A1a in the Supplementary Material (SM), is characterized by a dry summer and singular rainfall events during the winter, adding up to a mean annual precipitation (MAP) slightly less than 100 mm. We simplify the actual rainfall regime by considering periodic precipitation regimes with annual periodicity, + = P t P t ( 1) ( ), and an intra-annual square-wave pattern of the form where k is the number of rainfall events occurring at times t t , , k 1 … limited to the first half year, t t T is the equal duration of any rainfall event, A is the rainfall intensity and t ( ) θ is the unit step function [ x x x x ( ) 0 for 0, ( ) 1 for 0] θ θ = < = ≥ . An example for such a rainfall regime is given in Fig. A1b in the (SM). We vary the number of rainfall events (k) and their intensity (A) keeping the product kAT constant and equal to MAP 100[mm] = . This choice allows us to compare the responses of terraced riverbeds to precipitation regimes ranging from many weak rainfall events to a few strong events.
The response of a terraced riverbed to a particular rainfall regime is studied by calculating the year-end values of the biomass variables, χ = B B ( ) ij j i , at long times, from which the biomass distributions χ B ( ) j along the tradeoff axis χ for all terraces = … j N 1, , , are constructed. The model equations were solved numerically using forth order Runge-Kutta method. The number of species (m) was 500 and the system included 9 terraces. The daily precipitation was given as an input file according to Eq. (5). The simulations were integrated to sufficiently long times to ensure convergence to steady-state solutions. Numerical values for the model parameters used in these simulations are displayed in Table 1.

Data availability
There is no data available to this paper since it presents a theoretical work.