Highlights

  • Template and internal feedbacks interact to determine wetland distribution.

  • The relative strength of the two drivers varies with annual hydrological regime.

Introduction

A landscape is a mosaic of biotic and abiotic patches, and the configuration of this mosaic has potential consequences for ecological functioning (Strayer and others 2003; Dormann and others 2007; Ludwig and others 2007). Two important drivers produce this patchiness: (1) a preexisting physical template and (2) self-organization induced by internal feedbacks (Deblauwe and others 2008; Rietkerk and van de Koppel 2008; Liu and others 2014). A typical example of top-down control is physical constraint. For example, vegetation can be influenced by local topography through its effect on local microclimate or soil type, leading to template-dictated vegetation patchiness (Monger and Bestelmeyer 2006; Ropars and Boudreau 2012). Template-induced vegetation patchiness can also be generated by spatial distribution of resources (for example, water, nutrients, organic matter) (Kohler and others 2012; Manolaki and Papastergiadou 2012).

Biological patchiness can also arise from internal feedbacks in the absence of template heterogeneity. An example of patchiness arising from internal feedbacks is self-organization, which has been studied extensively in the past decade (for example, HilleRisLambers and others 2001; Rietkerk and others 2002; van de Koppel and others 2005; van de Koppel and others 2008; Weerman and others 2010). Regularly patterned vegetation (variously referred to as spots, bands, labyrinths, and gaps) occurs in many ecosystems around the world (Rietkerk and others 2002; Deblauwe and others 2008). The causal mechanism for regular spatial patterning often involves scale-dependent feedbacks. For example, resource concentration beneath vegetation improves local condition for plants, but results in resource depletion at distance, leading to local positive feedback and long-range negative feedback (Klausmeier 1999; Rietkerk and others 2002). Without the long-range negative feedback, local feedbacks alone can still generate spatial heterogeneity, for example scale-free patchiness, which can be described by a power-law distribution (Scanlon and others 2007; Kéfi and others 2007). Here, we define broadly the process of biological patch formation by internal feedbacks only (local positive feedback alone or both local positive feedback and long-range negative feedback) as ‘self-organization.’ ‘Self-organization’ investigated in this study refers to local positive feedbacks, without invoking long-range negative feedbacks.

Studies of self-organized patterns have focused on otherwise homogeneous landscapes, where the underlying template effect is absent or negligible; however, most real landscapes are heterogeneous and biological patchiness is a product of combined effects by both local feedbacks and preexisting template heterogeneity. Studies of spatial patterns created by interactions between the physical template and self-organization are still rare, however (Scheffer and Carpenter 2003). Sheffer and others (2013) studied a rock–soil mosaic system and found that when density of rocks was high, plant dispersion depended upon the distribution of rocks. Only when the open soil area between rocks was large enough, did plants form self-organized regular patterns. In addition, the relative importance of these drivers in organizing biological patchiness may shift with environmental conditions. For example, the importance of local facilitation, relative to competition, increases in response to increasing abiotic stress in forming subalpine and alpine plant communities (Callaway and others 2002). Recently, Dong and others (2017) have showed that the interaction between local feedbacks and template heterogeneity gave rise to nutrient spatial patterns in stream surface water across successional time in a stream ecosystem. Toward late succession as nutrients became more limiting, the local feedbacks played a progressively more important role in structuring spatial patterns of nutrient concentration.

Local facilitation and template heterogeneity also determine the responses of spatial patterns to environmental changes, including sudden discontinuous transitions of ecosystem state when environmental change reorganizes internal feedbacks that created the biological patchiness in the first place (Rietkerk and others 2004; Heffernan 2008). For example, mussel beds of ordered pattern are more resilient to wave disturbance than are random or clustered patterns (van de Koppel and others 2008) and spatial self-organization of mussel communities provides a refuge from predation in rocky intertidal systems (Robles and Desharnais 2002). Self-organized salt marsh ecosystems are resilient to strong wave attack on short timescales, but self-organization increases vulnerability of salt marsh to wave destruction on longer timescales (van de Koppel and others 2005). Meanwhile, theoretical studies suggest that underlying physical heterogeneity interferes with internal feedback-driven trajectories of ecosystem state change (van Nes and Scheffer 2005). How template heterogeneity interacts with concurrent self-organization to affect ecosystem resilience is still poorly understood. In one important study, Bonachela and others (2015) imposed an extra layer of soil water spatial heterogeneity induced by termite mounds onto a feedback model that described the self-organization of aboveground biomass and found that the savanna ecosystem became more resilient to a changing environment with preexisting spatial heterogeneity of soil moisture. Thus, it is critical to investigate how the heterogeneous template and internal feedbacks simultaneously shape biological patchiness in real landscape, which will enable us to better predict the response of ecosystems to environmental perturbations.

Desert streams are highly variable in both space and time (Fisher and others 1982; Grimm and Fisher 1989; Stanley and others 1997), making them ideal systems to investigate the dynamics of underlying drivers for biological patchiness. In particular, water permanence (defined as the proportion of time in a year a given site in the stream is covered by surface water) varies greatly in space (Stanley and others 1997), and the spatial distribution of riverine wetlands along desert streams is highly sensitive to the water permanence (Heffernan and others 2008; Dong and others 2016). However, this template may not solely determine the distribution of wetlands. Local positive feedbacks between vegetation and channel stability during floods affect vegetation persistence (Phillips 1990; Heffernan 2008), and the sediment deposition in wetland patches may enhance their resilience to drying (Heffernan and others 2008). The relative importance of the two drivers—the existing spatial heterogeneity of water permanence and local positive feedbacks—in determining wetland spatial patterns is yet to be evaluated.

Here, we addressed three questions: (Q1) What are the relative contributions of local feedbacks and the hydrogeomorphic template in shaping the spatial distribution of riverine wetlands in a desert stream? (Q2) How does the relative contribution vary with hydrological conditions? Finally, (Q3) how do template heterogeneity and local feedbacks influence wetland abundance and distribution in Sycamore Creek? To investigate these questions, we used a lattice simulation model parameterized by data on macrophyte distribution and abundance collected in a 12-km section of Sycamore Creek, Maricopa County, Arizona, USA. The data were collected annually during a 9-year period (2009–2017), encompassing the wide range of hydrological variability that characterizes this region. The model is phenomenological and does not represent individual feedbacks via distinct model terms (processes). Instead, ecological mechanisms underlying inter-annual vegetation dynamics are inferred from inter-annual variation in best-fit model parameters and their correlation with the annual hydrological regime.

Methods

Site Description

The study was carried out in Sycamore Creek (Figure S1), a tributary of the Verde River located 32 km northeast of Phoenix in the northern Sonoran Desert of central Arizona, USA. The stream drains a catchment of 505 km2. Mean annual precipitation ranges from 39 in lowlands to 51 cm year−1 in the mountains and exhibits winter and summer maxima with high inter-annual variability. The stream dries frequently, with surface water absent in many parts of the channel especially in summer, but water permanence varies substantially in space (Stanley and others 1997). Winter floods scour away most of the biota and initiate ecosystem succession in the following spring (Fisher and others 1982). The system experienced a state change around year 2000, with bare sand and cobble substrates colonized by algae transforming to a state with abundant riverine wetlands (Heffernan 2008), concentrated in a narrower valleys with more permanent surface flow (Heffernan and others 2008). Since then, wetlands (stream reaches dominated by vascular hydrophytes) have developed extensively in the system and their abundance and spatial distribution vary substantially among years in response to inter-annual variation in the hydrological regime (Dong and others 2016).

Data Collection

Data were collected from the 12-km mainstem of Sycamore Creek between 600 m and 700 m elevation each June from 2009 to 2017. The hydrology of the 9 years can be described in four categories (Figure 1; Table S1): very wet, wet, dry, and very dry. 2010 was a very wet year, which experienced a 100-year flood with a peak discharge of 439 m3 s−1 on January 22 and another large flood in late March 2010 (68 m3 s−1). All wetlands were removed by those events (that is, no visible aboveground plant parts). The years 2009, 2013, and 2017 were wet, with two or three medium-sized floods (peak discharge of 30–220 m3 s−1; Figure 1A). The years 2011, 2014, 2015, and 2016 were dry, with peak discharge of 10–70 m3 s−1 (Figure 1A). The driest year in the study period was 2012, with only one small flood (2.4 m3 s−1) in the mid-December. Such high hydrological variability is typical for Sycamore Creek and for other streams of the American Southwest.

Figure 1
figure 1

(A) Cumulative mean daily discharge (from October 1 of the previous year to the survey date of the next year) of the 9-year study period (2009 to 2017). Each abrupt increase in discharge represents a flood. (B) Relationship between cumulative mean daily discharge in flood days and wetland cover in 9 years and transitions among years.

Full size image

Each year on June 15 (± 3 days, near the time of peak biomass), the stream was surveyed from upstream to downstream and the location and size of each wetland patch was recorded. We used the five most abundant wetland species as indicator wetland species: Equisetum laevigatum, Paspalum distichum, Schoenoplectus americanus, Juncus torreyi, and Typha domingensis. We defined a wetland patch as a patch with at least one of the five indicator species present in an area greater than 4 m2, and we did not measure/record patches smaller than this size. If the gap between two patches was less than 1 m, the two patches were treated as one continuous patch. We used handheld GPS devices (with a resolution of 5 m) to record the upstream location of each patch, measured patch width and length with a tape measure, and recorded the dominant species within each patch.

Water permanence data were collected by E. Stanley (University of Wisconsin, Madison, WI, USA; Stanley and others 1997). She and her colleagues walked the same 12-km stretch of the stream every month for a consecutive 22 months between 1988 and 1989 and recorded the spatial extent of surface water. We used that dataset to calculate the proportion of the time in the 22-month period when surface water was present for each location (resolution is 1 m) along the stream. (Therefore, the water permanence is from 0 to 1.) In the model, we applied the same water permanence data collected in 1988–1989 to each year in our study period. We assumed that water permanence is a relatively stable feature of the geomorphic template, as it is controlled by drainage area and proximity of bedrock to the sediment surface (Figure S2).

Model Description

The model represents the stream states that are influenced by the hydrogeomorphic template (water permanence along the channel) and local feedbacks on transitions between gravel and vegetated state. The model is phenomenological and does not represent individual feedbacks via distinct model terms (processes); however, inter-annual variation in model parameters and their correlation with hydrological regime provides a basis for inferring underlying mechanisms influencing vegetation transitions. We used this model to fit wetland patch distribution for each year to assess the relative importance of template and self-organization. We used inter-annual variation in model parameters and its relation to annual hydrological regime to infer the mechanisms of the template and feedback effects. Finally, we manipulated the parameters for template and for feedback to explore the effect of these drivers on wetland abundance in the system.

The model was simulated on a one-dimensional lattice. Although we collected two-dimensional patch data (both length along the stream and width perpendicular to the channel), here we use only length data to simplify the model. The total number of 1-m-long cells was 12,000. Each cell could be in one of the two states: vegetated or empty (gravel, no plants). In the same study system, Heffernan (2008) demonstrated that a vegetated patch could be considered an alternative stable state to a gravel patch; that study had a spatial grain of about 20 m of stream length. Consistent with this scale of self-organization, we do not consider individual cells to exhibit alternative stable state dynamics. Instead, the local reinforcing interactions at the individual cell scale give rise to the bistability at the patch scale. We used water permanence (wp) as the hydrogeomorphic template, as it was found to be the major variable determining wetland abundance and distribution in the study system by previous research (Dong and others 2016).

In this model, the transition probability for a cell to shift from gravel to vegetated state, w01, is determined by two events: local dispersal and successful establishment. The latter event is further modified by the local water permanence condition (wpi, the water permanence value in cell i):

$$ w_{01} = q\varepsilon_{{\text{yr}}} {{wp}}_{i}^{{\beta_{{\text{yr}}} }} $$
(1)

where q is the percentage of cells that are vegetated in the neighborhood of the empty cell. These vegetated neighboring cells provide propagules and/or seeds to disperse to the empty cell. The neighborhood size was assumed to be 10 m (10 cells), with 5 m upstream and 5 m downstream. We tested for asymmetric influence from upstream and from downstream neighboring cells and found that it did not influence the results of this study. The key reasons might be that dispersal is mainly through vegetative reproduction, instead of seed dispersal, which will otherwise rely on flow. Although subsurface flow is also unidirectional, it is very slow (< 1 m h−1), so its effect on the direction of vegetative expansion is limited. Additionally, the wetland species we examined mostly grow along the bank of the stream, instead of directly in flowing surface water. We also tested a range of neighborhood sizes (from 2 to 40 m) and found the results to be robust with respect to all sizes tested. We assumed that all dispersal occurred locally because our test of a model including global dispersal did not improve model performance. The rate of propagule establishment or seed germination, ε, is the probability that a propagule dispersing from a neighboring vegetated cell successfully lands on an empty cell or a seed successfully germinates when the empty cell has water permanence of 100%. We assumed that propagule dispersal and establishment occur throughout the growing season; accordingly, water permanence represents the average hydrological condition of the growing season. After the propagule is established (or a seed germinates), the further growth and survival of a plant are affected by the hydrological condition of the local site. We divided the plant colonization into two stages: (1) germination (seedling establishment) ε and (2) plant growth and survival from seedling. We assumed that hydrology has different effects on these two stages. The first stage—germination (seedling establishment)—often occurs in late January/early February, which coincides with the time of winter floods. Seedling emergence has been shown to be significantly negatively influenced by sediment addition and flooding (Peterson and Baldwin 2004). The successfully established seedlings enter the second stage—plant growth and survival, which coincides with an extended period of declining stream flow. In this stage, success probability is reduced if the site has lower water permanence and dries out quickly, an assumption supported by previous study of the system (for example, Dong and others 2016). To describe this effect, we used a power-law formulation, with parameter β: The larger the reduction in plant establishment by local environmental condition, the greater the value of β.

The transition probability of a cell in the lattice changing from vegetated to gravel state, w10, is described by two variables: the mortality rate of an individual plant growing by itself without other plants nearby (that is, not in a wetland patch) and the reduction in the mortality rate caused by growing within a plant patch. Mortality rate is expressed as

$$ w_{10} = m_{\text{yr}} \left( {1 - q_{i}^{{{{f}}_{\text{yr}} }} } \right) $$
(2)

where myr is the year-specific mortality rate in the absence of local facilitation and q is the proportion of neighborhood cells that are vegetated. We modeled the effect of local facilitation using power-law formulation, taking advantage of the flexibility of the shape of the function. Parameter fyr (fyr > 0) is indicative of the year-specific unit strength of local facilitation: An fyr value closer to zero suggests stronger unit local facilitation. \( q_{i}^{{{{f}}_{\text{yr}} }} \) is indicative of realized mean local facilitation by taking into consideration the actual amount of wetland plants in a given year. Mechanisms that could produce such facilitation include positive effects of plant biomass (density) on channel stability during floods events (Heffernan 2008) or/and on water-holding capacity and soil water content in the rooting zone during dry periods (Caldwell and others 1998). Our multi-year study allows us to evaluate the likely importance of these mechanisms. If local facilitation is realized via the mechanism of density-dependent channel stabilization during flooding, we expect the effect to be more pronounced in a wetter year with more frequent floods. If facilitation is caused by increased water-holding capacity in the plant rooting zone, we predict a stronger effect in drier years.

We used periodic boundary conditions to avoid edge effects. (We tested the potential cascading effects from a boundary by creating scenarios in which there was an upper and/or lower boundary that was fixed in the bare or vegetated state. We did not detect boundary effect on model results, likely explained by the spatial heterogeneity of the template and the local spatial dispersal in the model compared to the size of the system.) The lattice was updated asynchronously in order to approximate a continuous-time process (Durrett and Levin 1994). The model randomly selected 500 cells as vegetated cells for the initial condition. We tested the effect of number of random cells (ranging from 200 cells to 11,500 cells) in the initial condition on model performance and found that they converged to a statistically identical state (Figure S3). During each simulation, cells were randomly selected for update, with N individual cell updates at one time step, until a statistical steady state was reached, that is, the total number of vegetated cells stabilized. (On average, it takes ~ 300,000 iterations to stabilize; in each iteration, the model updates the state of 12,000 cells.) We used three metrics, namely patch size distribution, gap size distribution, and match of cell state (gravel and vegetated) between modeled and observed result, to quantify goodness-of-fit of the model and determine the best parameter set for each year. To process the observed data for each year, we lined up the wetland patch distribution by position along the 12-km stream, divided evenly into 12,000 cells. Each cell was labeled as 1 (with plant) or 0 (without plant), according to the wetland distribution. We calculated the observed patch size and gap size distribution for each year to compare with model results. Each error was estimated by the mean squared deviation between data and predicted values normalized by the data variance. We also compared the observed and modeled results at each cell across the entire 12-km stream and computed the percent of mismatches in cell state (presence/absence of plants). The goodness-of-fit model was evaluated as the mean of error of patch size distribution, error of gap size distribution, and the percentage of mismatches. We explored the parameter space within reasonable ranges and calculated the total error for each parameter combination. The combination with minimum total error was selected as the best-fit parameter set (Table S2). Using this procedure, we determined the best-fit parameter set for each year.

With this model structure, there are two inherent system-wide equilibria: (1) when the whole landscape is vegetated, probability of transition to an unvegetated state w10 (Eq. 2) equals 0 everywhere and (2) when the whole landscape is gravel (unvegetated) w01 (Eq. 1) equals 0 everywhere. These trapping states occur when m is very low or very high, but the model can reach intermediate equilibria in wetland extent (a steady state with a mix of gravel patches and wetland patches) with other parameter values. In Appendix A, we present a detailed description of the model using the standard ODD (overview, design concepts and details) protocol (Grimm and others 2006).

Analyses and Numerical Experiments

We calculated year-specific cumulative mean daily discharge during flood days (the days when the mean daily discharge > 1 m3 s−1) for 2009–2017 to describe the hydrological state of Sycamore Creek for each year. The conclusions derived in this study are not sensitive to the formulation of this hydrological metric (Figure S4). In Sycamore Creek, the ‘wetness’ of a year is highly correlated with the frequency and magnitude of winter floods; a year with frequent and/or large winter floods is very likely to have more extensive and longer-lasting surface flow, and thus be considered a wet year. We analyzed the statistical correlation between parameter values (mortality, m; seedling establishment rate, ε; template effect parameter, β; and mean realized local facilitation corrected for the effect of total cover, \( \tilde{f} \)) and the annual cumulative flood flows. The year-specific realized local facilitation is calculated by \( \frac{{\mathop \sum \nolimits_{i = 6}^{11995} q_{i} f_{\text{yr}} }}{11990} \) (starting with i = 6 because a neighbor size of 5 cells is used in the model). The value of this index is influenced by total wetland cover. To correct for the compounding effect of total cover on the relationship between realized local facilitation and annual cumulative flood flow, we further calculated the value of local facilitation under the given total cover when plants are randomly distributed (zero spatial autocorrelation among vegetated cells). We then analyzed the statistical correlation between cumulative flood flow and the realized local facilitation corrected for the effect of total cover (that is, \( \frac{{\mathop \sum \nolimits_{i = 6}^{11995} q_{i} f_{\text{yr}} }}{11990} \) under actual distribution minus the value under random distribution). The analysis of correlation was carried out for the cumulative flood discharge on both the original scale and on the log scale. (Log transformation was applied to reduce the effect of the single extreme wet year of 2010 and also to make the values in dry years more dispersed.)

We used a numerical procedure to explore the effects of local facilitation and template heterogeneity (that is, water permanence) on wetland cover among years varying in hydrological conditions (see ODD in Appendix A). To determine the effect of template heterogeneity on wetland distribution and abundance (to address Q3), we constructed three templates: (1) actual template: using the actual water permanence; (2) random template: created by randomizing the spatial distribution of the actual water permanence to achieve zero spatial autocorrelation; and (3) homogeneous template: using the mean of the actual water permanence as water permanence for all 12,000 cells. We compared the distribution of wetlands with these different templates under each year’s hydrological regime. To infer the effect of local facilitation on wetland distribution and abundance, we explored the response of wetland extent to variation in the value of fyr. To do so, we used a range of values that center on the best-fit value for that year. We then ran the model to a statistical steady state to determine its equilibrium state of wetland occupancy.

Results

Model Performance and Patch Size Distribution

The abundance of wetlands along the 12-km mainstem of Sycamore Creek varied greatly in the 9-year study period, ranging from less than 1500 m in 2010 to about 5000 m in 2012 (Figure 1B; Figure S5). In general, drier years featured higher wetland cover, but the wetland abundance could vary widely under very similar hydrological conditions (Figure 1). On average over 9 years, 23% (varies between 18 and 33%) of the 12,000 cells switched their state between vegetated and gravel states from 1 year to the next (Table S3).

The mean error across 9 years was 0.18 (Table 1; averaging over the study period, 29% of the mean error was contributed by mismatch in gap size distribution, 18% by mismatch in patch size distribution, and the rest 53% from mismatch in cell state—gravel and vegetated state). For all 9 years, the tail of the patch size distribution followed a power law (R2 > 0.9; Figure 2) (that is, values greater than some minimum, xmin; Clauset and others 2009). The value of the exponent ranged between 1.20 and 1.82 in different years. (A larger exponent means a faster drop in tail, indicating a higher proportion of smaller patches.) Other than the two outlier years 2013 and 2014 (the 2 years following the driest year), the exponent of power-law distribution of wetland patches was significantly linearly correlated with the total wetland cover of each year (Figure 3).

Table 1 Year-Specific Goodness-of-Fit of the Model and Parameter Values in the Best-Fit Models
Full size table
Figure 2
figure 2

Goodness-of-fit of the model for patch size (upper two rows) and gap size distribution (lower two rows) between 2009 and 2017. The observed patch size distributions (the solid black dots) were fitted with power-law distributions (the solid gray line), with γ being the exponent. a on the x-axis is a given patch size (or gap size), and y-axis is the number of patches (or gaps) greater than the patch size a (or gap size a) in the system. Gray open circles are results from best-fit model for each year.

Full size image
Figure 3
figure 3

The exponent of power-law distribution of wetland patch size (see Figure 2) was significantly linearly correlated with the total wetland cover of each year, except for 2013 and 2014—the 2 years following the severe drought periods (2011 and 2012). Dashed lines indicate the 95% confidence interval of the linear regression model.

Full size image

Relationship Between Local Facilitation, Template Heterogeneity, and Hydrology

Both mortality m and germination rate ε decreased with discharge (Figure 4A, B); that is, wetter years were associated with lower mortality and lower rate of successful seedling establishment. However, neither relationship was statistically significant (p = 0.13 for mortality and p = 0.29 for seedling establishment rate; Figure 4A, B and Table S4). Template effect (β) increased with discharge, suggesting stronger template effect on plant establishment rate in wetter and more frequently flooded years (p = 0.01 on the original scale and p = 0.1 on the log scale of annual cumulative flood flow; Figure 4C and Table S4). The correlation between f value and cumulative flood flow was not statistically significant (p = 0.80). The strength of realized mean local facilitation (after correction for the effect of total wetland cover) increased with discharge values, and such trend was statistically significant (p = 0.034 on the original scale and p = 0.021 on the log scale of cumulative flood flow; Figure 4D and Table S4).

Figure 4
figure 4

Relationship between the estimated strength of different processes (indicated by year-specific best-fit parameter values: m—mortality rate; ε—rate of success germination of seeds or establishment of propagules; β—effect of template on plant survival; and qf—realized mean effect of local facilitation on plant survival corrected for the effect of total wetland cover) and the annual accumulative mean daily discharge in flood days (mean daily discharge > 1 m3 s−1) between 2009 and 2017.

Full size image

The form of the template effect varied among years (Figure 5A). In drier years, the relationship between template effect on plant establishment rate was more concave, suggesting that the establishment rate was most sensitive to changes in water permanence in places where surface flow does not last long (lower water permanence). In contrast, in wetter years the curve was less concave, indicating approximately equal level of sensitivity of establishment rate across the water permanence range, and more widespread occurrence of drying stress. Averaging 1-\( {\text{wp}}_{i}^{{\beta_{\text{yr}} }} \) (as defined in Eq. 1) over the entire stream can be used to quantify the mean template effect. For the period between 2009 and 2017, the mean template effect (quantified by 1-\( {\text{wp}}_{i}^{{\beta_{\text{yr}} }} \)) is 0.063, 0.25, 0.024, 0.024, 0.072, 0.070, 0.070, 0.11, and 0.070, respectively. Averaging over 9 years, the grand mean is 0.084—that is, water permanence on average reduced the plant survival rate by about one-tenth across the entire stream.

Figure 5
figure 5

(A) Effect of water permanence gradient on the rate of successful propagule establishment between 2009 and 2017, that is, wp βi yr. (B) Year-specific mortality rate as a function of percent patch cover in the neighboring cells (upstream 5 m and downstream 5 m) in each year, that is, \( 1 - q_{\text{yr}}^{{f}} \). The different years are color-labeled in four categories: blue—very wet year (2010); red—very dry year (2012); green— wet years (2009, 2013, and 2017); and magenta—dry years (2011, 2014, 2015, and 2016). (Color figure online)

Full size image

The form of the effect of local facilitation on mortality rate also varied among years, but much less markedly than the template effect (Figure 5B). All f values were greater than 1, with mean of 2.21, describing a strong nonlinearity between effect of local facilitation and patch size. To reduce mortality rate by 50% through the mechanism of local facilitation requires about 75% plant cover in neighboring cells (Figure 5B). The mean realized strength of local facilitation is estimated by \( q_{i}^{{{{f}}_{\text{yr}} }} \) (as defined in Eq. 2; without correcting for the effect of total wetland cover), and the values for the study period between 2009 and 2017 were 0.18, 0.074, 0.31, 0.30, 0.23, 0.14, 0.13, 0.20, and 0.14, respectively. The mean over all 9 years was 0.19—that is, local facilitation reduced the mortality rate by about one-fifth on average over the study period across the entire stream. As a result, the mean effect of local facilitation on plant mortality was about twice as much as the mean effect of template heterogeneity on plant colonization rate.

Effect of Local Facilitation and Template Heterogeneity on Ecosystem Changes

Effect of Local Facilitation

When unit strength of local facilitation was weakened in our numerical experiments by increasing f value, the system eventually shifted to the gravel state (Figure 6). On the other hand, sufficiently strong unit local facilitation led to a system to be completely covered by wetland plants (Figure 6). An intermediate level of unit local facilitation sustained a heterogeneous landscape with a mix of gravel and vegetated patches (Figure 6). Wetland cover in years with stronger unit local facilitation (smaller f values) was more sensitive to changes in the strength of local facilitation (Figure 6J).

Figure 6
figure 6

Effect of local facilitation (greater f, weaker facilitation) on the steady-state wetland cover across the 12-km mainstem of Sycamore Creek (AZ) between 2009 and 2017 (A–I). (J) relationship between Δf (defined as the difference in f value between the maximum f when the system is entirely covered by wetlands and the minimum f when the system transitioned to the gravel state, that is, the shaded windows in A–I) and best-fit f value for each year.

Full size image

Effect of Template Heterogeneity

The presence of template heterogeneity was much more important than the distribution, whether real or random, for explaining wetland extent. In the homogeneous template, the system eventually collapsed to a gravel state regardless of hydrological conditions (Figure 7). Except for year 2010, there was no significant difference in wetland cover between a system imposed with a random template and the one with an actual heterogeneous template (Figure 7)—that is, the degree of spatial autocorrelation did not influence wetland cover. In 2010—the year with strongest template effect (Figure 4C)—wetland cover was significantly higher when the underlying template was random compared to the case of actual spatial heterogeneity (Figure 7B).

Figure 7
figure 7

Effect of physical template (water permanence gradient along the stream), that is, homogeneous template (‘H’), the actual template (‘A’), and random template (‘R’), on steady-state wetland cover in Sycamore Creek (AZ) between 2009 and 2017.

Full size image

Discussion

Abrupt shifts in ecosystem state are rooted in internal feedbacks that also favor creation of the biological patchiness (Rietkerk and others 2002; Heffernan 2008). Meanwhile, underlying physical heterogeneity can interplay with internal feedback-driven trajectories of ecosystem state change (van Nes and Scheffer 2005). Our 9-year study demonstrated that internal feedbacks and template heterogeneity operated simultaneously to influence wetland distribution and that the relative importance of each varied greatly among years depending upon hydrology. The average effect of local facilitation during the 9-year study period was about twice that of the mean template effect and about 11 times greater in the driest year of 2012.

Inter-annual Variability of Wetland Cover

The correlations between plant mortality m and seedling establishment rate ε and the year-specific discharge were not statistically significant, implying weak relationship of hydrological conditions and wetland abundance. This is contradictory to the previous findings, which suggested the paramount effect of hydrological regime on wetlands in arid streams (for example, Casanova and Brock 2000; Dong and others 2016). Such contradiction is likely caused by the lumped nature of the model parameters: Wetter years will favor higher flood-caused mortality and lower drought-caused mortality, and vice versa in dry years and same for seedling establishment rate ε. Overall, in drier years, seedling establishment is widespread (Figure 4B) and less sensitive to the water permanence (Figure 4C); but established plants die more often (Figure 4A). In wet years, seedling establishment is less common (Figure 4B) and depends strongly on the gradient of surface water permanence (Figure 4C); but once established, plants die rarely (Figure 4A). These suggest that high flows may inhibit seedling emergence (Peterson and Baldwin 2004), likely because atmospheric oxygen levels are required for germination and early seedling growth (Karssen and Hilhorst 1993). The cumulative effects of these tend to restrict the abundance of wetlands overall in wet, more flooded years (Figure 1B).

In years of similar discharge regime, total wetland cover in the system still varied by more than 100% (Figure 1). Such a large inter-annual variation is likely related to the legacy effect of hydrology of the proceeding year(s) affecting seed banks and propagules, a key mechanism that enables persistence through variable hydrological conditions in desert streams and rivers (Capon 2003). Here we assume that the occurrence of wetland plants is not limited by seed banks in the system; instead, the successful expression of seed banks—consequently, the distribution of wetland plants—is limited by hydrology conditions, as supported by studies carried out in other dryland streams (for example, Bagstad and others 2005). Several lines of evidence from our study implied that the major drought period of 2011–2012 may have caused a shift in the state of system that required several years to return to the pre-drought state. Immediately following 2012, the year 2013 had a disproportionally higher wetland cover compared to 2009 and 2017, whose hydrology was very similar to that in 2013 (Figure 1). In two consecutive post-drought years (2013 and 2014), the power-law distribution describing the patch size distribution exhibited a much larger exponent than predicted by total wetland cover (Figure 3), indicating a higher proportion of relatively small patches than predicted by the wetland cover.

A high proportion of small patches could be a consequence of weaker local facilitation (Table 1) and/or the stronger role played by local factors (for example, sediment texture, soil types, reach types). Sycamore Creek is frequently flooded, and in such hydrologically flashy systems, the plant community is predominantly shaped by a landscape-scale flood regime (Capon 2005). In less flooded years, dry periods stimulate germination and modify oxygen availability in the substrate, leading to the general trend of higher plant cover in drier years (Figures 1B, 3). The over-expression of wetlands in 2013 following the two dry years was likely a legacy effect through stimulation of seed banks. In these times, the negative biotic factors (for example, competition) become more important, which diminish the local positive facilitation assumed in the model, resulting in a higher proportion of small patches (Figure 3). Lacking specific information on performance of individual plant species, however, inference of competition remains conjecture. In contrast, the system bounced back with little legacy effect after the 100-year flood in 2010 (Figures 1B, 3). This is consistent with many other studies showing that drought—as opposed to floods—is more likely to trigger a shift in system state (for example, Capon 2005; Bogan and Lytle 2011).

Local Facilitation

The mechanism of local facilitation takes diverse forms in harsh environments. In intertidal zones, plants reduce flow erosion and enhance sediment stabilization (Bruno 2000). In subalpine and alpine plant communities at high elevations, vegetation maintains local temperatures and protects from strong winds (Callaway and others 2002). In arid and semiarid terrestrial ecosystems, plants increase local water and/or nutrient availability, forming ‘island of fertility’ (Schlesinger and Pilmanis 1998). Studies of plant communities show that local feedbacks vary in response to abiotic stress (for example, Callaway and others 2002). In wetlands in arid drainages, our study suggests stronger overall local facilitation in drier years (Figure 4D). We expected that local facilitation could be realized (1) if plants modify soil texture via root zone action that facilitates their survival in dry periods and/or (2) if channel stabilization by plant roots facilitates resistance to flood scouring during frequently flooded, wet periods. The model result suggested that in Sycamore Creek, the mechanism to resist drying seems to be more effective and influential than the mechanism to resist flood disturbance (Figure 4D). This inference is derived from the result of stronger realized local facilitation in drier years than in wetter years (Figure 4D). However, the inference is confounded by the effect of inter-annual variation in patch size on the strength of realized local facilitation—wetter years were characterized by smaller patch sizes (Figure 2), lowering local facilitation. Additionally, we did not find statistically significant correlation between unit strength of local facilitation and annual hydrological condition. In the very dry years (2011 and 2012), we observed considerable plant mortality in our annual survey, suggesting that facilitation by modifying local temperature and soil moisture may not be effective enough against severe droughts, a conclusion shown elsewhere (for example, Casper 1996). These lines of evidence are indicative of an overall stronger role of local facilitation in increasing resistance of the wetland community to hydrological disturbances, and they do not distinguish the dominant mechanism of local facilitation in maintaining wetland plants under extreme hydrological regimes.

Vegetation establishment is strongly related to broader-scale hydrogeomorphic feedbacks, not just local processes (Larsen 2019). Studies of riverine wetlands have shown reciprocal adjustments between fluvial landforms and vegetation dynamics (Corenblit and others 2007). For example, during floods in vegetated streams, plants often control morphogenetic processes in channels and on floodplains by promoting stabilization and sediment accretion (May and Gresswell 2003). This sediment stabilization modifies hydrological disturbance patterns along fluvial corridors, driving possible positive feedback on plant survival (Bendix and Hupp 2000). The combined effects of microscale soil–vegetation positive feedback and large-scale vegetation–climate positive feedback have been shown to amplify nonlinearity, generating alternative stable states (Janssen and others 2008). Broader-scale feedbacks occur at timescales greater than 1 year, which requires a legacy effect reflected in the model to account for nonlinear dynamics (Zweig and Kitchens 2009). Our model did not account for potential legacy effect by precedent year(s)’ hydrology (for example, on the propagule vitality), which could confound our interpretation of parameter values: They could have integrated the effect of hydrology from both the current year and precedent year(s).

Template Heterogeneity

The template effect is stronger in organizing the plant community during more flooded years (Figure 4C). In desert streams, sites of higher water permanence often occur where an open valley pinches into a canyon (Stanley and others 1997; Dent and others 2001). Conversely, sediment deposition occurs as water velocity is reduced at more open sections of the stream (Grant and Swanson 1995). Hence, the spatial gradient of water permanence may correlate with the erosion–deposition gradient. If so, sections of higher water permanence would be associated with stronger erosion and lower water permanence with greater sediment deposition. In wetter, more flooded years, the spatial variability of erosion–deposition and that of the water permanence would synchronize to produce an even larger difference in plant survival rates among sites of different valley floor widths, resulting in the greater template effect. However, our finding of a significant correlation between the template effect and hydrology relies heavily on the result of one year, 2010—a historically wet year that shows much stronger template effect than any other years. Additionally, the correlation between template effect and hydrology was statistically significant (p = 0.0069) on the original scale of the annual flood flow, but only marginally significant (p = 0.098) on the log scale of the annual flood flow. Log transformation reduced the influence by the ‘outlier’ data point, that is, the extreme wet year of 2010 (Table S4), rendering the correlation insignificant. A more robust analysis of multi-year lags will require a longer time series to capture diverse sequences of annual hydrological conditions.

We demonstrated the essential role of template heterogeneity in sustaining wetland cover: Under our model, in a homogeneous setting, the system eventually shifts to gravel state regardless of hydrological condition (Figure 7). In reality, vegetation–flow–sediment feedbacks could link vegetation abundance to drying stress in ways that would counteract this theoretical trajectory. In a heterogeneous environment, however, locations with higher water permanence could serve both as a ‘refuge’ and as a ‘reservoir’ in that they allow plants in these locations to persist when faced with high environmental stress and to recover when the stress is reduced. In our study, the variable controlling spatial heterogeneity was maintained at a constant (same mean water permanence in the ‘actual,’ ‘random,’ and ‘homogeneous’ scenarios). Therefore, the difference in ecosystem response to environmental change is solely caused by template heterogeneity. Similarly, by fixing mean soil moisture, Bonachela and others (2015) imposed a layer of soil–water heterogeneity induced by termite mounds onto a feedback model that described the self-organization of aboveground vegetation. They found that spatial heterogeneity of soil moisture increased the resilience of savanna ecosystems to climate change.

Spatial autocorrelation generated differences in wetland cover only in the wettest year of 2010, when the wetland cover with the random template was twice that of the actual heterogeneous template (Figure 7B). The random template (zero spatial autocorrelation) is significantly patchier than the actual template (Figure S2). Such patchiness in water permanence could be only effectively harnessed for plant establishment when both the strength of local facilitation and the template effect are high, a condition satisfied only in wettest years (for example, 2010). It is worth noting that in reality, in addition to spatial heterogeneity of water permanence, streams are notoriously heterogeneous: characterized by heterogeneity of channel substrata, a longitudinal riffle–run–pool arrangement, and variation in nutrient availability. Many of these factors may affect vegetation distribution. Spatial heterogeneity in multiple dimensions could make the ecosystem response to environmental perturbations much more complex, potentially explaining why most direct evidence for ecosystem bistability has been only from models and lab experiments (Schröder and others 2005).

Averaging over the 9-year study period, the effect of local facilitation was about twice that of the template effect in Sycamore Creek, and about 12 times in the very dry 2011 and 2012 (Figure 5). The relative strength of local facilitation and template heterogeneity affects the response of ecosystems to environmental changes. Positive feedbacks amplify small deviations, which can destabilize the system globally (Kéfi and others 2016). Template heterogeneity also alters ecosystem resilience, the strength of which depends on the dispersion among patches (van Nes and Scheffer 2005). Local facilitation studied here can be interpreted as a type of dispersion, through which the consequences of local hydrological conditions are propagated to influence others nearby, resulting in heterogeneous systems behaving like single homogeneously mixed systems. When dispersion is strong and template heterogeneity is spatially autocorrelated, a domino effect occurs in the system’s response to environmental perturbations (van Nes and Scheffer 2005). To conclude, our study shows template effect and internal feedbacks operate simultaneously to influence the distribution and abundance of wetland plants in a dryland stream. The relative effect of each varied from year to year, influenced by hydrology. Template and self-organization together interacted with a complex suite of other factors to generate residual year-to-year variation, with an implied strong role for legacy effects. Desert stream ecosystems and their wetland communities are likely to be more sensitive to drought, predicted to become more frequent and severe in this region due to global climate change, than to flooding.