Coupled local facilitation and global hydrologic inhibition drive landscape geometry in a patterned peatland

Introduction Conclusions References


Introduction
The structure and function of natural ecosystems are shaped by complex interactions between biotic and abiotic processes acting at different spatial scales.In resource-limited environments, these interactions can give rise to self-organized, patterned landscapes (Rietkerk and Van de Koppel, 2008;Dyskin, 2007).Archetypal examples exist in arid and semiarid ecosystems (Foti and Ramirez, 2013;Saco et al., 2007;Scanlon et al., 2007;Klausmeier, 1999;Mabbutt and Fanning, 1987) and peatlands (Prance and Schaller, 1982;Eppinga et al., 2009;Larsen and Harvey, 2011).These patterns range from regular mosaics with characteristic length scales to scale-free patterns exhibiting heavy-tailed patch size distributions (von Hardenberg et al., 2010).While both types of patterns typically signify the resource-limited nature of their respective environments, the primary biotic and/or abiotic processes that dictate the evolution of regular and scalefree landscapes are thought to be considerably different.Regardless of their driving mechanisms, patterned landscapes create ecological heterogeneity and thus help maintain biological diversity (Kolasa and Rollo, 1991) and productivity, and increase system resiliency (van de Koppel and Rietkerk, 2004).Given their reliance on a critical resource (e.g., water, nutrients, or both), the presence of self-organized patterning also suggests that even subtle disturbances or environmental changes can lead to catastrophic shifts in ecosystem states (Kefi et al., 2007(Kefi et al., , 2011;;Reitkerk et al., 2004).Therefore, understanding the mechanisms that govern the development and maintenance of landscape pattern is crucially important = 0.98) based on the results of numerical simulation of surface water flow using a hydrodynamic model (SWIFT2D; Schaffranek, 2004).
to conserving their ecological attributes, particularly as the exogenous drivers are disrupted by climate change and largescale anthropogenic landscape modification.
The ridge-slough landscape in the Everglades (Florida, USA) is a patterned landscape in which two distinct vegetation communities, ridges and sloughs, comprise a selforganized landscape mosaic (Fig. 1a).Ridge patches occupy higher soil elevations ( z ∼ 25 cm in the best-conserved landscapes; Watts et al., 2010) and are dominated by the emergent sedge Cladium jamaicense, while sloughs contain a mix of submerged, floating leaved, and emergent plants.One of the most striking features of the landscape is its clear anisotropic spatial structure, with ridge patches elongated in the direction of historic surface water flow (SCT, 2003;Larsen et al., 2007).Hydrologic modifications have led to the loss of this patterning over large areas of the historical ridge-slough landscape since the turn of the 20th century (Watts et al., 2010;Nungesser, 2011).Rapid conversion to randomly oriented, isotropic ridges, and even towards monotypic sawgrass or cattail (Typha domingensis) landscapes, is attributed mainly to compartmentalization and attendant hydrologic modification, along with agricultural nutrient enrichment (Wu et al., 2006;Light and Dineen, 1994;Gaiser et al., 2005).
Adverse ecological impacts associated with reduced slough extent and connectivity (SCT, 2003) have led to pattern maintenance and restoration being focal points for restoration planning and assessment.A prerequisite for successful ecological restoration efforts is a clear understanding of pattern genesis mechanisms and the time frame for that process to occur.Despite multiple plausible mechanisms, these feedbacks remain poorly understood (Cohen et al., 2011) principally because of difficulties associated with observation and experimental manipulation at large enough spatial and temporal scales, a constraint that focuses attention on models of pattern genesis and degradation.
Several hypotheses have been proposed to explain the ridge-slough patterning and have produced models that are capable of generating elongated patches.In all cases, pattern evolution in modeled landscapes responds to water flow (and flow direction), but it remains unclear which flow attributes and/or processes govern the process.Moreover, it has been shown elsewhere that different combinations of processes may generate similar patterns (Eppinga et al., 2009), making inference of a single dominant mechanism challenging.Models of ridge-slough pattern genesis have invoked differential sediment transport (Larsen et al., 2011;Lago et al., 2010), nutrient redistribution (Ross et al., 2006), and biased subsurface flow induced by anisotropic hydraulic conductivity (Cheng et al., 2011) as driving mechanisms.Using a process-based numerical model, Ross et al. (2006) showed that differential evapotranspiration rates in higherelevation soils may lead to the concentration of dissolved nutrients (particularly phosphorus), suggesting that this nutrient redistribution, alone or in combination with sediment and nutrient transport, may generate the ridge-slough pattern in the Everglades.According to Larsen andHarvey (2010, 2011) and Lago et al. (2010), ridge-sloughlike flow-parallel patterns develop as the heterogeneous flow regimes caused by vegetation resistance and local elevation differences recursively dictate differential peat accretion, sediment transport, and erosion.Cheng et al. ( 2011) incorporated anisotropic hydraulic conductivity in a scaledependent feedback (advection-diffusion) model to demonstrate the evolution of flow-parallel, elongated vegetation bands.Crucially, however, all of these models generate pattern by coupling local positive feedbacks with a distal negative feedback at some intermediate distance imposed either explicitly (Cheng et al., 2011) or implicitly (e.g., erosion and sediment transport by flow; Larsen andHarvey, 2010, 2011;Lago et al., 2010).
As we discuss below, although these models yield patches that are elongated in the direction of flow, the modeled landscapes lack several other critical pattern metrics unique to the conserved portion of the ridge-slough system.With the exception of Larsen andHarvey (2010, 2011), all other models generate patterns that are strikingly regular and with a characteristic separation distance (e.g., Cheng et al., 2011), a property that appears to be absent in the Everglades ridgeslough mosaic.Furthermore, our recent analyses (Casey et al., 2015) of the ridge-slough landscape have shown that the patches in this landscape are strongly scale-independent, suggesting that the negative feedback that stabilizes patch expansion in these landscapes is a global, rather than a distal, scale-dependent mechanism.
An alternate explanation for ridge-slough patterning that includes a global negative feedback is the self-organizingcanal (SOC) hypothesis (Cohen et al., 2011;Heffernan et al., 2013).In this conceptualization, spatially anisotropic patterning emerges from global constraints on patch (both ridges and sloughs) expansion created by feedbacks between landscape geometry, water flow, and hydroperiod, which controls peat accretion (and therefore affects landscape geometry).In short, the SOC hypothesis proposes that elongated patches develop as the landscape incrementally adjusts its spatial geometry (ridge density, size, and shape) to optimize the discharge competence (i.e., ability to convey water; Heffernan et al., 2013;Kaplan et al., 2012;Cohen et al., 2011), and that this pattern may evolve without sediment or nutrient redistribution.Decrease in discharge competence primarily results from high ridge density but can be further intensified by reduced patch elongation (which lowers the probability of longitudinally connected sloughs).Both scenarios yield a global increase in hydroperiod for a given boundary flow (Kaplan et al., 2012).Increased hydroperiod, in turn, makes conditions more favorable for ridge-to-slough transitions, which decreases ridge density, lowering hydroperiod and ultimately tuning landscape pattern to hydrology.
The core assumption in the SOC hypothesis is that patch elongation occurs because changes in patch density affect discharge competence anisotropically.Specifically, expansion of ridge patches parallel to flow has a neutral effect on discharge competence (and thus hydroperiod), while patch expansion orthogonal to flow has a strong negative effect on discharge competence.Kaplan et al. (2012) demonstrated this feedback with a hydrodynamic model of surface water flow across randomly generated patterned landscapes of varying anisotropy, finding that hydroperiod decreased exponentially with increasing patch anisotropy.Using a patchscale analytical model, Heffernan et al. (2013) demonstrated strong feedbacks between the soil elevation at any given location and an adjacent location perpendicular to flow, but no such feedback parallel to flow, lending support to the central mechanism of the SOC hypothesis.However, that analytical model was limited to two patches, where flow in one cell was directly controlled by flow in the only other cell.Further testing the SOC hypothesis requires evaluating the potential for this mechanism to generate anisotropic patterning at the landscape scale.
In this study we implemented the local and global feedbacks described in the SOC hypothesis (Cohen et al., 2011;Heffernan et al., 2013) in a stochastic cellular automata framework to model temporal evolution of the ridgeslough pattern.Transition probabilities between ridge and slough states were driven by hydroperiod (a global negative feedback) and local facilitation.Both isotropic and anisotropic neighborhood kernels were implemented, and simulations were performed under different combinations of local-facilitation strength and degree of anisotropy to investigate the role of each process in pattern development.Simulated ridge-slough landscapes were then compared to a suite of statistical and geostatistical characteristics that characterize the ridge-slough patterning observed in the bestconserved remnants of the Everglades.

Hydrologic model
We first expanded the hydrodynamic modeling procedure outlined by Kaplan et al. (2012) to calculate hydroperiods for landscapes over a wide range of ridge densities (%R) and anisotropy values (e).Steady-state flows were simulated for 960 synthetically derived ridge-slough landscapes, with %R and e ranging from 10 to 90 % and 1.0 to 6.0, respectively, using a spatially distributed numerical flow model (SWIFT2D) (Schaffranek, 2004).For each landscape, rating curves were created by applying a series of constant head boundary conditions (BCs) at the up-and downstream model domain boundaries, assuming uniform flow and no-flow BCs at the domain lateral boundaries.A 20-year (1992-2012) daily time series of modeled flows (see Kaplan et al., 2012) was then used along with landscape-specific rating curves to calculate daily surface water elevations and resulting hydroperiod in each landscape.Finally, a polynomial function response surface fitted to modeled data yielded a meta-model of hydroperiod (HP, the fraction of time a location is inundated) as a function of %R and e (Fig. 1b).

Cellular automata model
Our cellular automata (CA) model consists of two states: ridges and sloughs.While the natural system contains variations in these states (e.g., wet prairie communities that can persist in short hydroperiod sloughs; Zweig and Kitchens, 2008), they are likely transitional states between sloughs and ridges and were not included.Tree islands were likewise neglected in the CA model.While tree islands are critically important to the Everglades landscape, they represent only approximately 3 % of the total landscape area, and their emergence and maintenance is thought to be controlled by different mechanisms than those explored here (Ross et al., 2006;Wetzel et al., 2009).
System states in the ridge-slough landscape are differentiated by two primary characteristics: vegetation and soil elevation (Watts et al., 2010).Our model assumes that, when a cell transitions from one state to the other, vegetation and soil elevation attributes are updated immediately.The probabilities that govern transitions between states are dictated both by a global feedback based on hydroperiod and a localfacilitation effect of neighboring cells.A schematic of the model framework shows the recursive algorithm that generates landscape pattern (Fig. 2).

Local facilitation
Local facilitation of patch expansion was modeled based on the similarity of neighboring cell states to the cell state under transition.That is, the probability of a cell changing state is locally controlled by the neighborhood of adjacent cells.Ecologically, these effects may arise due to plant propagation (vegetative and reproductive), changes in primary production at patch edges that change peat accretion rates, and potential abiotic factors, such as nutrient and sediment transport mediated by flow (Ross et al., 2006;Larsen et al., 2011).Local facilitation, λ, ranges from 0 to 1 and decreases exponentially with distance (d), such that where j = 1, 2, . . ., n (number of neighbor cells); x is the state of each neighbor cell (x = 1 if the same state as the center cell, 0 otherwise); and k is an exponential decay parameter.Large values of k in Eq. ( 1) indicate that the localfacilitation effect decreases sharply within a short distance, indicating that the immediate neighbor cells contribute most of the local facilitation.Small k values, in contrast, connote a larger area of effect with all neighbors contributing more uniformly to the local facilitation.A circular neighborhood was employed with the radius determined such that the cumulative distribution function of Eq. ( 1) for all surrounding neighbor cells within that radius exceeds 99 % (e.g., Foti et al., 2012;Scanlon et al., 2007).Note that Eq. ( 1) assumes an isotropic neighborhood effect; i.e., neighboring cells influence adjacent cells regardless of direction.Directional bias on local facilitation, which may occur in response to the direction of flow, was also imposed by varying the local-facilitation decay rate, k, as a function of the angle between the center cell and its neighboring cells: where α j is the absolute azimuth angle between the longitudinal direction and a line connecting the cell and its neighbor (expressed as 0 ≤ α ≤ 90 for each quadrant of the neighborhood kernel), k x is the exponential decay parameter in the east-west (E-W) direction, and k y is the exponential decay parameter in north-south (N-S) direction.In Eq. ( 2), the ratio k x : k y describes the directional bias of the local-facilitation effect.

Transition probabilities
The HP meta-model is the hydrologic foundation of the CA model (Fig. 2), with HP variation creating the global negative feedback that drives changes in ridge-slough configuration.
Landscape patterns conducive to efficient drainage (lower %R, higher e) lower hydroperiod, increasing the probability of a slough pixel transitioning to ridge, while landscape patterns that inhibit drainage (higher %R, lower e) decrease that transition probability.Crucially, this hydroperiod effect is manifest at the domain scale (i.e., globally).We assume landscape HP as the global driver instead of a local parameter because of the extremely low relief (ca.0.003 %) of the Everglades which creates highly uniform water levels over the entire wetland system.Transition probabilities between ridge and slough were modeled as the linear combination of local effects (i.e., by surrounding neighbors) and global effects (i.e., controlled by landscape HP), expressed as where λ R and λ S are the local facilitation for ridge and slough, respectively, and HP t = 0.87 is the target hydroperiod, based on the range of expected HP for well-conserved landscapes in the Everglades (Givnish et al., 2008;McVoy et al., 2011;Cohen et al., 2011).
The transition probability formulations (Eqs.3 and 4) are identical to those used by Foti et al. (2012), with one key difference: these authors imposed the global negative feedback mechanism by directly setting a target vegetation density, whereas here a target hydroperiod is implemented.This formulation is based on observations of the temporal dynamics of change in the ridge-slough landscape, which suggest that ridge density can change quickly towards a landscape that is dominated by either ridge or slough based on hydroperiod (Nungesser, 2011).Furthermore, this construction allows for variable ridge density driven by HP, which, in turn, is controlled by the density and anisotropy of patches (Kaplan et al., 2012).Setting a target HP therefore explicitly considers bidirectional interactions between hydrology and landscape geometry, allowing for future modeling based on perturbations to hydrological forcing.

Model domain and parameterization
Simulations begin with a randomly generated, 3.5 km × 3.5 km landscape composed of 10m cells with low %R (5-15 %), following the suggestion that ridges formed out of a slough matrix (Bernhardt and Willard, 2009).At each time step, transition probabilities for each cell are calculated based on Eq. ( 1) through ( 4).This probability matrix is used to determine transitions between ridge and slough cells, producing a new landscape in the following time step.Based on the new landscape configuration (i.e., different %R and e), a new HP is calculated (Fig. 1), yielding a new global feedback function based on Eqs. ( 3) and ( 4).Landscape configuration thus changes iteratively (Fig. 2), eventually reaching a dynamic equilibrium when %R, e, and HP have stabilized.
The local feedback mechanism is dictated by the magnitudes of k x and k y as well as their ratio.Increasing magnitudes of k x and k y indicate greater "spatial immediacy" of neighbors (Scanlon et al., 2007) in the respective directions, while the ratio k x : k y describes the magnitude of anisotropy in local facilitation.We explored effects of magnitude and ratio of these two parameters by simulating landscapes using seven different values of k y from 0.10 to 0.40 (corresponding to a circular neighborhood with radius from ca. 20 to 70 m) and seven levels of anisotropy in local facilitation (k x : k y ratios from 1.0 to 4.0).With increasing k x : k y ratio, the shape of the local-facilitation effect becomes more elliptical with elongation in the N-S direction as local feedback in the E-W decays sharply with increasing distance.Simulations were replicated a minimum of six times for each combination of k y and k x : k y values to ensure that equilibrium landscape characteristics were controlled by model parameters rather than random initial conditions.We note from preliminary simulations that, for a given k y , a maximum value of k x : k y existed, beyond which any increase in k x severely restricted patch expansion in the E-W direction.This resulted in extremely narrow ridges (1-2 cells wide) that drove e values towards ∞, leading to model instability because of the feedback of e on hydroperiod.For example, for k y = 0.2, steady state could not be attained for k x : k y > 3.5, while for k y = 0.35 steady state could be achieved only for k x : k y ≤ 2.0.Therefore, for each k y magnitude the final simulations were performed only for the k x : k y ratios that consistently resulted in a stable equilibrium.

Landscape pattern metrics
Simulated ridge-slough landscapes from the CA model were compared, both qualitatively and quantitatively, to observed ridge and slough patterns in the best-conserved portion of the Everglades (referred hereafter as "reference landscapes") (McVoy et al., 2011;Watts et al., 2010).Thirteen reference landscapes from a study of Everglades landscapes by Nungesser (2011) were augmented by eight additional refer-ence landscapes presented by Casey et al. (2015).All pattern metrics were based on analyses of rasters created with a 10 m cell resolution from vector maps (Rutchey et al., 2005).
Comparisons between simulated and reference landscapes were based on seven statistical and geostatistical characteristics: overall %R, e, correlation length scales (semivariograms) parallel and orthogonal to historical flow, distribution of patch sizes, perimeter area fractal dimension (PAFRAC), and landscape characteristic wavelength (periodicity).
Patch density was calculated as ridge area divided by total domain area.Patch anisotropy was estimated as the ratio of the major (parallel to flow) and minor (orthogonal to flow) ranges of indicator semivariograms (Deutsch and Journel, 1998); the correlation length scales inferred from these spatial ranges were also of interest to ensure that the model predicted realistic patch geometry.Distributions of patch sizes, which follow power-law scaling in the reference landscapes, were evaluated for goodness of fit in comparison to the Pareto (power law) distribution using Monte Carlo tests (Clauset et al., 2009).The fractal dimension (PAFRAC), which measures patch shape complexity, was calculated from the fitted slope between patch area and perimeter.Finally, patch periodicity was evaluated using radial spectrum (r-spectrum) analysis, which extracts the spectral components of the landscape pattern as a function of possible wave numbers (i.e., spatial frequency divided by domain size).R-spectra, which are used to identify the characteristic wavelength and directional components of regular patterns, were obtained using two-dimensional Fourier transforms following methods outlined by Couteron and Lejeune (2001).
Our primary motivation in comparing model results to reference landscapes for these metrics is to elucidate the nature of local interactions required to create pattern that is statistically consistent with the best-conserved portions of the extant Everglades (i.e., elongated, flow-parallel ridge patches).We therefore applied a multi-criteria objective function to quantify agreement between simulated and reference landscapes for the pattern metrics listed above.Simulated landscapes received a score of 1 for each of the seven metrics that fell within the range of values observed in reference landscapes, and the sum was used as an integrated measure of pattern agreement (IMPA; maximum IMPA = 7.0).A mean IMPA score was calculated for each parameter set (i.e., the average IMPA score of 6 simulation replications), allowing us to identify combinations of parameters that yielded patterning concordant with the reference landscapes.While the IMPA approach is useful for identifying model parameters that generate patterning consistent with the well-conserved ridge-slough mosaic, it is also important to note that some pattern metrics have been found to be ubiquitous across the present-day Everglades.For instance, power-law scaling and lack of a characteristic separation distance (aperiodicity) are observed in all parts of the ridge and slough system -regardless of their state of degradation (Casey et al., 2015).This suggests that a reliable ridge-slough landscape model should produce landscapes with these criteria for all possible parameter combinations (i.e., both isotropic and anisotropic model formulations), whether the simulated landscapes resemble the conserved or degraded ridge-slough mosaic.

Results
Hydrodynamic modeling of the discharge competence and landscape hydroperiod suggested that %R exerts dominant control on the landscape hydroperiod, while anisotropy is of secondary importance.Consequently, the ridge density of simulated ridge-slough landscapes showed greater sensitivity to changing hydroperiod regimes than the patch anisotropy.Patches in simulated landscapes with symmetric local facilitation (k x = k y ) were always isotropic (e < 1; Fig. 3), with low ridge density (ca.33 %) compared to reference landscapes.Implementing directional bias (k x : k y > 1.0) in local facilitation resulted in anisotropic landscapes with higher %R.However not all k x : k y configurations yielded patterning geostatistically similar to the reference landscapes.For example, for all landscapes where k y = 0.1 (effective neighborhood radius = ca.70 m), regardless of the k x : k y ratio, patches were highly diffuse, extending across the entire domain, with spatial correlation lengths both parallel and orthogonal to flow that were much larger than in the reference landscapes (Fig. 3).As k y increased, patches became more distinct and aggregated, with patch geometry that better resembled the reference landscapes both qualitatively and quantitatively.
Simulated landscapes for k y = 0.2 (effective neighborhood radius = ca.40 m) and k x : k y of 1.0, 2.0, 2.5, and 3.5 show a strong qualitative resemblance to an example reference landscape (Fig. 4a), particularly at higher values of  k x : k y (≥ 2.5).High values of k x : k y also resulted in reference landscape similarity for semivariogram ranges (Fig. 4b) and pattern anisotropy, with e values within the range observed in reference landscapes (2.5-5.0)(Nungesser, 2011;Kaplan et al., 2012).
The distribution of ridge areas in the simulated landscapes showed significant support for power-law scaling following Clauset et al. (2009) (i.e., we cannot reject the hypothesis that the distribution differs significantly from power law at the 0.1 level).Power-law scaling was consistent across model realizations, regardless of %R and e (i.e., for all parameter combinations, Fig. 5a), which agrees with the observed patch size distribution in both conserved and degraded ridge and slough landscapes (Casey et al., 2015).Moreover, the values of the fitted exponent (β) were strikingly similar between reference (β = 1.73 ± 0.09) and simulated landscapes (β = 1.73 ± 0.0.11), and this exponent did not vary significantly among landscapes simulated using with different k x : k y ratios -a result also echoed by the similarity of fit-  Casey et al. (2015).
The perimeter-area scaling of patches showed that the modeled landscapes were highly fractal (Fig. 5b), as evidenced by the linear relationship between log (perimeter) and log (area) (slope > 0.5; Foti et al., 2012).However, the scaling relationship for the reference landscapes deviate from the linear function, indicating their non-fractal nature.Although a linear-scaling relationship seemed to hold acceptably within a certain patch-size range (Fig. 5b), Casey et al. (2015) found that the perimeter-area scaling in real ridgeslough landscapes was better explained by a quadratic function over the entire size range (i.e., without any cutoffs).This indicates that the larger patches in real landscapes are increasingly more complex as opposed to the simple power functions followed by the modeled landscapes that suggest an equal degree of shape complexity for all patch sizes.Finally, similar to the r-spectra in reference landscapes, our simulated landscapes exhibit the clear absence of a characteristic wavelength (i.e., r-spectra maxima at nonzero distance) (Fig. 5c).This suggests that both the real and simulated landscapes are aperiodic (i.e., not regularly patterned) (Casey et al., 2015), a feature consistent with global rather than scaledependent patch constraints.
Summarizing the statistical and geostatistical properties of simulated landscapes (Fig. 6, symbols), in comparison with values observed in reference landscapes (shaded region), illustrates the relatively narrow parameter space over which model outputs match the conserved (i.e., elongated, N-S oriented) patterning.Both %R and e increase nonlinearly with increasing k y and k x : k y (Fig. 6a-b), and the slope relating e and k x : k y is highest when k y is large.The E-W semivariogram range (Fig. 6c) declines exponentially with k x : k y for all k y , while the N-S range (Fig. 6d) remains relatively flat.In combination, our results suggest that only a subset of simulated landscapes meet the multiple conditions observed in reference landscapes (i.e., have IMPA = 6), while power-law scaling, patch complexity, and aperiodic patterning were evident in all simulations, regardless of parameterization.Notably, when k y is > 0.30 (effective neighborhood radius ≤ ca. 25 m), most (88 %) of the simulated landscapes aligned with reference landscapes (IMPA = 6.0), even when k x : k y is relatively small (i.e., 1.5 < k x : k y ≤ 2.5).In contrast, with lower k y (e.g., k y = 0.15, equivalent to a neighborhood radius of ca.45 m), only 24 % of simulations met the multi-criteria objective, and they required the highest value of k x : k y to do so.

Testing the self-organizing-canal hypothesis
A clear understanding of the processes underlying development of ecological patterns is integral to all ecosystem man- agement and restoration.In the Everglades, venue for one of the largest and most ambitious ecosystem restoration efforts in history, the specific focus on landscape pattern as a restoration objective underscores the urgency of the process-pattern link.Identifying the suite of necessary and sufficient processes to create and maintain pattern will aid in prioritizing hydrologic restoration goals.Although multiple hypotheses exist for explaining the ridge-slough pattern, most of them attribute the development of these landscapes to one dominant process.The self-organizing-canal hypothesis (Cohen et al., 2011;Heffernan et al., 2013), on the other hand, ascribes pattern formation and maintenance to reciprocal feedbacks between landscape pattern and hydrology.Moreover, evidence of a strong feedback between pattern and hydroperiod (Kaplan et al., 2012) lends support for the SOC.Primacy of this mechanism vis-à-vis nutrient enrichment or sediment redistribution -and we note here that these mechanisms are not mutually exclusive -would imply markedly different water management objectives, specifically emphasizing flow volume sufficient to ensure appropriate hydroperiod vs. water level management or creation of episodic high velocity The dominant spatial feature in the ridge-slough landscape is patch orientation with flow.As a minimum criterion, models that fail to produce flow-oriented elongation are clearly insufficient explanations for pattern development.Our results suggest that the SOC mechanism can create patterning consistent with the best-conserved ridge-slough landscape, but only when local facilitation is directionally biased in the direction of flow.This suggests the SOC alone is an insufficient mechanism.Previous work comparing static landscapes (Kaplan et al., 2012) showed that anisotropy exerted strong control on hydroperiod, but our model results suggest that in a dynamic landscape changes in patch density (%R) occur more quickly than changes in patch shape (e).As such, %R provides the dominant control on HP (Fig. 1b), while anisotropy plays a secondary role.When local feedbacks are isotropic, global negative effect on the landscape pattern due to hydroperiod selects for landscapes with %R at ca. 30 %, far below the value in the best-conserved pattern, and is unable to generate patch anisotropy.
Recently, Heffernan et al. ( 2013) used an analytical model to explore the SOC, demonstrating that ridge and slough elevation divergence occurs spontaneously at some discharge levels and that the impact of a given cell on adjacent cells orthogonal to flow is far larger than parallel to flow.In short, pattern arises solely due to feedbacks between hydroperiod and discharge competence (i.e., capacity to convey water), which is controlled by the configuration orthogonal to flow.To reconcile these findings with our model results, we note that water flow in the Heffernan et al. (2013) model is limited to two flow paths, where occlusion of flow in one cell (e.g., due to peat accretion there) must, of necessity, force water through the other.In contrast, our model comprises a relatively large domain with hundreds to thousands of possible flow paths, weakening the influence of flow occlusion in any given cell on global hydroperiod.As a result, the role of anisotropy on discharge competence is diminished, and ridge density impacts on hydroperiod dominate.
We also note that the HP in this study is estimated assuming steady-state flow conditions (as calculated over the 20-year period of record in Kaplan et al., 2012) and does not represent the possible effects of temporal fluctuations in flow that occur in the Everglades ecosystem.In order to test whether a fluctuating hydrological regime would drive elongated ridge formation under isotropic local facilitation, a variable hydrology scenario was also implemented in the CA model based on reported variation in mean flow into Lake Okeechobee over a 65-70-year cycle (Enfield et al., 2001).However, the simulations driven by cyclically varying hydrology coupled with isotropic local facilitation did not drive ridge elongation in the resulting landscapes (i.e., e = 1) and yielded low values (< 30 %) of %R due to the recurring high-HP events.These initial simulations suggested that variation in HP affected %R much more than e and was not sufficient to drive anisotropic patch evolution.

Multi-metric model performance
As noted above, the model developed here does create compelling ridge-slough patterning given anisotropic localfacilitation effects.The plausibility of multiple model mechanisms, including those presented here, for creating floworiented elongation suggests that additional landscape characteristics are necessary for evaluating model performance.The additional proposed pattern metrics (ridge density, anisotropy, autocorrelation range, patch size distribution, fractal dimension, and periodicity) provide a more nuanced and comprehensive basis on which to compare model outputs to real landscapes.This approach is similar to Larsen and Harvey (2010) wherein multi-metric comparisons between modeled and real landscapes were made, but includes new potentially relevant pattern metrics.While it was beyond the scope of the current work to compare the multiple existing models of the ridge-slough landscape, we note that our model outputs agree reasonably well with observations in the best-conserved ridge-slough landscapes for all of the proposed metrics.
Among the most important differences between our model and others for the ridge-slough pattern is invocation of an inhibitory feedback that operates globally rather than at a characteristic spatial scale.Constraints to patch expansion are induced at the entire domain scale, and not over local and/or intermediate scales, as is the case in Ross et al. (2006), Lago et al. (2010), Cheng et al. (2011) and Larsen and Harvey (2011).The principal reason for invoking a global rather than intermediate feedback is the inherent difference between focusing on hydroperiod/water depths, which are reasonably uniform over large areas, and flow velocity or solute redistribution, which are more spatially heterogeneous.Global feedbacks have been widely invoked to understand and simulate vegetation patterning (e.g., Scanlon et al., 2007;Foti et al., 2012), and induce three pattern features that merit particular attention: power-law scaling of patch areas, high fractal dimension (i.e., highly crenulated patch edges), and the absence of a characteristic pattern wavelength (which would be expected in regular patterning).Our simulations closely matched observed power-law scaling of patch areas (including the scaling parameter, β) and, perhaps most importantly, the absence of a characteristic pattern wavelength implying no regular landscape periodicity.These landscape properties are exhibited in both well-conserved and degraded ridgeslough landscapes in the Everglades (Casey et al., 2015).That our model consistently reproduces them suggests that global inhibition is integral to ridge-slough pattern evolution.However, that it does so across all model runs, even those that clearly fail to reproduce credible patterning, means that these metrics are necessary -but not sufficient -for discriminating the anisotropic pattern genesis processes.
Pattern geometry -including flow-oriented elongation, which is the sentinel feature of this landscape -is strongly controlled by the local-facilitation function in our model.The remaining three metrics (%R, e, and semivariogram ranges) were used to conclude that only a subset of parameterizations for inducing local feedbacks yielded landscapes with geostatistical properties in agreement with reference landscapes.For example, %R and e of modeled landscapes fall well outside the reference values at low k y or small k x : k y ratios (Fig. 6a, b).Likewise, modeled landscape semivariogram ranges tend to be well above what is observed in the real systems (Fig. 6c, d) for these parameterizations.While the particular mechanisms that induce anisotropic facilitation are, as yet, unclear (see below), we can at least conclude that some parameterizations of local and global controls can satisfy all diagnostic metrics (Fig. 6e).The spatial range of the extant pattern is controlled in the model by k x and k y , which control the distance over which local facilitation acts in the E-W and N-S directions, respectively.While statistically compelling landscapes can be simulated using several parameter values, a synthesis of model performance (Fig. 6e) suggests all matching landscapes have k y > 0.15, corresponding to a local-facilitation kernel that extends, at maximum, 40 m in the N-S direction; this is further constrained in the E-W direction due to anisotropic facilitation (i.e., k x > k y ).Overall, of the 37 simulations, 22 % comported with all six metrics observed in the conserved landscapes, all of which required k x > k y and k y > 0.15.
Power-law scaling of patch sizes has been associated with vegetation self-organization in many landscapes (e.g., Scanlon et al., 2007;Kefi et al., 2007Kefi et al., , 2011) ) but has only recently been evaluated for the Everglades (Foti et al., 2012) and specifically for ridge-slough patterned landscape (Casey et al., 2015).Most studies of the ridge-slough landscape have emphasized the perceived regular nature of the pattern, including invocation of a pattern wavelength of ca. 150 m (SCT, 2001;Larsen et al., 2007;Watts et al., 2010).Notably, power-law scaling of patch area is incompatible with regular patterning because the basis of such patterning is the presence of distal negative feedbacks that truncate patch expansion at some particular spatial range (van de Koppel and Crain, 2007).It is therefore critically important that our simulated landscapes, wherein inhibitory feedbacks are global and not scale-dependent, exhibit this power-law scaling behavior (Fig. 5a).It is also notable that the landscapes follow such scaling regardless of ridge density or local-facilitation parameters, which is also observed in all reference landscapes with various patch densities (Casey et al., 2015).These results are consistent with the concept of robust criticality in ecological systems, where local spatial interactions lead to power-law clustering of patches well below the percolation threshold (Kefi et al., 2011;Vandermeer et al., 2008).While the generality of power-law scaling in both simulated and real landscapes limits this metrics utility as a model di-agnostic, it lends strong support for the primacy of scale-free processes underlying ridge-slough pattern formation.
Interestingly, patch complexity in the real ridgeslough landscapes revealed a non-fractal nature (nonlinear perimeter-area scaling; Casey et al., 2015).Since the simulated landscapes show a highly linear perimeter-area scaling and hence highly fractal patterns, this highlights one of the attributes of the ridge-slough landscapes that our model is not able to entirely reproduce.In contrast, Foti et al. (2012) recently reported that sawgrass patches, the dominant vegetation of the ridge-slough landscape in the Everglades, were fractals.However, their landscapes were analyzed at significantly coarser scale (40 m pixel) as opposed to 10 m pixel as in this study, which is likely to miss the finer scale crenulations in patch edges that increase the complexity.
The ridge-slough landscape is often described as exhibiting a repeating geostatistical pattern with a wavelength of 50-400 m in the direction orthogonal to flow (Larsen and Harvey, 2010;Lago et al., 2010;Cheng et al., 2010;Cohen et al., 2011).Spatial periodicity in patterned ecosystems has been attributed to the interplay between positive and negative feedbacks acting at different spatial scales.Short-range facilitation causes vegetation aggregation in dense clusters, but patch expansion is inhibited by some intermediate-range negative force acting at a specific distance.In this way, vegetation self-organizes into a periodic configuration (Rietkerk and Van de Koppel, 2008;von Hardenberg et al., 2010).Accordingly, the models presented by Ross et al. (2006) and Cheng et al. (2010) generate highly uniform, elongated patches that possess a clear periodicity.Surprisingly, however, the observed ridge and slough landscape appears to lack periodic spatial structure (Fig. 5c; Casey et al., 2015), suggesting there is no characteristic wavelength to the landscape.Even more surprising is that this aperiodic behavior is retained across a wide gradient of hydrologic modification.It is therefore notable that the model presented here lacks periodic spatial structure.The lack of landscape periodicity argues strongly against invocation of intermediate-range negative feedbacks.The observed pattern is more consistent with a global negative feedback that inhibits patch expansion across the entire landscape.
We note that the spatial scale (i.e., spatial resolution and extent) can strongly influence various landscape pattern metrics (e.g., Wu et al., 2002;Levin, 1992) we have used in this study.Geostatistical methods (e.g., semivariogram) are inherently affected by cell size (Lausch et al., 2013;Atkinson and Tate, 2000;Atkinson, 1993), while cellular automata models are also influenced by cell and neighborhood size (Pan et al., 2010;Ménard and Marceau, 2005;Chen and Mynett, 2003).Our modeling results and interpretations are based on 10 m grid size.While the minimum mapping unit (MMU) varies from 20 to 50 m (Nungesser, 2011;Rutchey et al., 2005), smaller features (< 10 m) are apparent in these mapping products.Setting raster and model resolution at 10 m captured the majority of perceivable features without

S. Acharya et al.: Coupled local facilitation and global hydrologic inhibition drive landscape geometry
requiring untenable computation times.The neighborhood size in our model is controlled by local-facilitation parameters k x and k y , which highlights that different neighborhood sizes produce patterns with remarkably different spatial attributes, and only a few parameter combinations can produce the patterns that are highly consistent with the reference ridge-slough landscape.

Mechanisms of anisotropy
Mechanisms of local facilitation in patterned landscapes are generally attributable to more than one biotic/abiotic factor, which can be difficult to measure or determine at landscape scales (Cohen et al., 2011).Our model yields novel insights about the role of a generalized local-facilitation process and its spatial extent in emergent ridge-slough patterns (e.g., local-facilitation effects confined to 40 m parallel to flow and even less perpendicular to flow).However, while the model creates compelling pattern based on the combination of global inhibition and anisotropic local facilitation, the mechanisms that induce anisotropic facilitation remain unclear.Several potential mechanisms exist.Flow may enable directional seed dispersal (i.e., hydrochory;Nilsson et al., 2010) particularly for sawgrass.Crucially, however, despite prolific seed production, most sawgrass reproduction is vegetative (Miao et al., 1998).Local anisotropic facilitation may also arise from sediment entrainment and deposition (Larsen et al., 2007), though we note that the invoked flow velocity effects to date have focused on inhibitory feedbacks (i.e., constraints on patch expansion), not local-facilitation effects.However, if deposition occurs preferentially downstream (e.g., at the tails of ridges) rather than at ridge edges, the cumulative effect would be anisotropic facilitation.Another mechanism posits lower phosphorus uptake efficiency in ridges than in sloughs, leading to longer uptake lengths (sensu Newbold et al., 1981), and thus further downstream transport of available P in ridges.While this mechanism remains untested, it comports with observations of substantial P enrichment in ridge soils compared with adjacent sloughs (Ross et al., 2006;Cheng et al., 2011) and could yield a directional stimulatory effect on sawgrass primary production.
While determining the mechanism that controls localfacilitation effects is clearly critical for successfully protecting and restoring landscape pattern, our work suggests that processes driving ridge-slough pattern development and maintenance may be represented by a generalized localfacilitation function and a global inhibitory feedback, potentially signifying a unifying explanation of ridge-slough pattern development.The model results presented herein provide the first test of ridge-slough simulations against a suite of expanded landscape-scale statistical and geostatistical properties, several of which strongly support inference of a dominant role for global feedbacks between pattern and hydroperiod in structuring this sentinel landscape.

Figure 2 .
Figure 2. Schematic representation of steps in the cellular automata model of ridge and slough pattern development.The upper central panel is a third-order polynomial surface of hydroperiod (HP) vs. anisotropy (e) and ridge patch density (%R) (R 2= 0.98) based on numerical simulation of surface water flow using a hydrodynamic model (SWIFT2D;Schaffranek, 2004).

Figure 3 .
Figure 3. Simulated landscapes for various kx : kyratios for (a) ky = 0.1 and (b) ky = 0.20.Note the increase in distinctiveness of ridge and slough patches with increasing magnitude of ky.

Figure 6 .
Figure 6.Mean values of statistical and geostatistical metrics in simulated landscapes (symbols) relative to the ranges observed in reference landscapes (shaded regions): (a) ridge density (%R), (c) patch anisotropy (e), (c) semivariogram ranges in the E-W direction (perpendicular to flow), (d) semivariogram ranges in the N-S direction (parallel to flow), and (d) average (IMPA) scores for all kx : ky combinations.