A Stochastic Approach toModeling Tidal Creek Evolution: Exploring Environmental Influences on Creek Topologies Through Ensemble Predictions

A numerical tool integrating hydrodynamics, vegetation growth, stochastic selection, and Monte‐Carlo simulation is developed to predict channel network evolution in a newly inundated back‐ barrier wetland. Focusing on the formation and evolution of tidal channels, model simulations suggest a clear tendency towards specific creek topologies with high probability despite creek evolution being subject to random perturbations. These creek topologies are primarily determined by wetland topography while the influence of vegetation is secondary. Ensemble results of Monte‐Carlo simulations identify threshold points leading to bimodal separation of the emergent creek networks characterized by different complexity and efficiency. This, in turn, assists the recognition of necessary human interventions to encourage the formation of complex and far‐reaching creek networks that can promote exchange of materials between estuaries and salt marshes and, hence, maintain salt marsh ecosystem structure and function. Plain Language Summary A large part of the world's low‐lying areas is heavily managed to protect lives, assets, and landscapes from marine‐derived natural hazards. However, climate change is putting unprecedented pressure on these managed environments, forcing the adoption of “no active intervention” or “managed realignment” strategies in areas where “hold the line” options cannot be justified due to financial constraints. This leads to the very likely possibility that unmaintained coastal defenses will begin to fail, raising questions as to the potential impacts on low‐lying regions as well as the options available to mitigate impacts through targeted interventions. Of particular interest is the question of whether it is possible to facilitate the shift of managed coastal wetlands towards intertidal environments in which tidal creeks form and evolve to support salt marshes. Combining the use of a numerical model and multiple model runs, we find that although the evolution of tidal creeks is subject to random disturbances, creek networks are likely to form with a clear tendency towards specific types. Grouped results of multiple, random simulations are useful in identifying interventions that are necessary for developing far‐reaching tidal creek networks. Influence of vegetation is found to be secondary compared to the wetland topography in terms of influencing tidal creek evolution.


Introduction
Adaptation of coastal areas facing climate change is a global challenge. Apart from the fact that about 10% of the world population lives in these low-lying regions (McGranahan et al., 2007) and that much infrastructure in these lowlands is critical to the local and global economy (e.g., Hsiang et al., 2017;Muis et al., 2015), our coastline also embraces natural environments providing habitats for wildlife and supporting recreational activities for humans Mudd et al., 2009;Temmerman et al., 2013). Consequently, these areas are commonly managed and engineered to reduce damage, loss of life, and environmental degradation caused by natural hazards originating from the sea. However, sea-level rise and increased storminess due to global warming (e.g., Sayers et al., 2015) are putting unprecedented pressure on these managed systems, rendering some of them uneconomical and unsustainable to maintain, in which case we adopt either a "no active intervention" or a "managed realignment" strategy (Dale et al., 2018;French, 2006;Turner et al., 2007).
The aim of this research is to examine the potential for reduced complexity modeling for exploring the constraints of wetland topography and vegetation on the nature of tidal creek evolution which plays a key role in the establishment of salt marsh plants and the maintenance of salt marsh stability (Kearney & Fagherazzi, 2016;Symonds & Collins, 2007). Wetland topography includes depressions that simply join a creek network ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. or concentrate flows and cause localized erosion, leading to creek expansion. Vegetation acts to moderate erosion and creek expansion by retarding flow and enhancing the resistance of the wetland to erosion. Knowledge of such constraints can provide guidance for resource and habitat managers in situations where "no active intervention" and "managed realignment" are defined strategic priorities.

10.1029/2019GL085214
To simulate creek evolution, D'Alpaos et al. (2005) combined a friction-dominated approximation for hydrodynamics (Rinaldo et al., 1999) and a random selection (Kirkpatrick et al., 1983) for headward growth of tidal channels. However, despite the randomness associated with the model, the results were regarded as deterministic, and the analyses were based on individual simulations. In contrast, more sophisticated shallow water equations have been used in combination with the standard advection-diffusion equation for calculating sediment transport to simulate the formation and development of tidal channels (Belliard et al., 2015(Belliard et al., , 2016Temmerman et al., 2005Temmerman et al., , 2007. Despite being able to achieve mass conservation and produce deterministic results, compared to the simplified model mentioned above, the use of shallow water equations requires significantly increased computational effort, leaving long-term simulations over large domains unfeasible (Mariotti, 2018).
Within the context of considering vegetation in the model, the presence of vegetation in D' Alpaos et al., 2007 was described by linking vegetation to local topography through a parametric relationship. This concept was also adopted in Kirwan and Murray (2007) and Mariotti (2018). In contrast, vegetation development was simulated by Temmerman et al. (2007) based on a thorough plant growth model considering initial establishment, lateral expansion, logarithmic growth of stem density, and plant mortality due to flow stress and inundation.
This research couples the Poisson hydrodynamic model proposed by Rinaldo et al. (1999), which requires much less computational effort, with the comprehensive vegetation growth model proposed by Temmerman et al. (2007) to provide a foundation for the simulation of tidal creek evolution that is based on the framework introduced in D' Alpaos et al. (2005).

Method
The major interactions between the hydrodynamics, vegetation growth, and creek evolution in the model are presented in supporting information Figure S1. A key element of tidal creek evolution module is that the erosional direction of the creek network in each time step is selected in a probabilistic manner. This is done via a method presented by D' Alpaos et al. (2005) where deterministically calculated bed shear stresses are used to influence the above-mentioned probabilistic selection. Due to the probabilistic nature of the model, no two simulations will yield the same creek network. As each model run results in a different and equally plausible creek, concentrating on a single output could be deceiving, potentially masking highimpact and low-probability creek formations. Hence, Monte-Carlo simulation (Metropolis & Ulam, 1949) is applied to obtain a statistical interpretation of the principles of tidal creek evolution. In the Monte-Carlo simulations, the models were run for 100 times each, with settings and external conditions of the model kept identical for each model run.
The model is applied to the RSPB Minsmere Nature Reserve which is located on the coast of Eastern England and fronted by a line of coastal sand dunes and sand and gravel beach ridges ( Figure S2). This site is selected because the management options for the sections of coastline are "managed realignment" or "no active intervention" (Table S1). The grid resolution is 10 m × 10 m throughout the domain. The topography of the model is defined based on data downloaded from EDINA DIGIMAP for the offshore areas and LiDAR data with 2 m resolution for the inland regions. Topographical features in the backbarrier wetland, such as existing drainage channels and embankments, are introduced in the model through manual adjustment.
Model runs are based on two types of initial site configuration: one with artificially created preexisting channels (dimensions given in Figure 1) within the site connected to a short breach channel at the Minsmere sluice (hereafter Case 1) and another with just the breach at the sluice (hereafter Case 2). The first condition is representative of a "managed realignment" intervention aimed at the production of a tidal creek network; the second one represents a "no active intervention" scenario with eventual barrier breaching. Cases 3 and 4,   Figure S2). Initial establishment of vegetation in the domain is modelled stochastically by multiplying an initial stem density of 10 NP/m 2 (number of plants per m 2 ) and a matrix of random numbers following a uniform distribution. Critical bed shear stress is set to 0.1 N/m 2 for grid cells with stem density less than 230 NP/m 2 , representing that of silt-very fine sand substrate which is a very common mix in lowland areas (Walker et al., 2004) and of areas undergoing early stage of vegetation establishment (Temmerman et al., 2007). The critical bed shear stress is increased to 0.3 N/m 2 beyond the threshold vegetation density of 230 NP/m 2 to reflect the soil reinforcement impact of vegetation on critical bed shear stress (e.g., Clark & Wynn, 2007;Dietrich et al., 1993). Values of the parameters of the model are detailed in Table S2. Geometry efficiency proposed by Marani et al. (2003) and complexity (ratio between the total channelized stream length and the radius of the smallest circle that covers the entire channelized area) are calculated to characterize the predicted creek networks. Detailed discussion around these two metrics can be found in Text S1. The developed model was applied to Hesketh marsh to examine the model's effectiveness in simulating a tidal creek network resulting from a recent managed realignment scheme in the Ribble Estuary, UK (see Text S2).

Statistics Reveal Bimodal Distribution of Creek Networks
Figure 1 includes eight creek networks randomly selected from the 100 networks created by Case 1, showcasing the diversity of the pattern of the predicted creek networks. Figure 2 shows the geometry efficiency and complexity of the 100 creek networks created by Case 1. The mean of geometry efficiency is 2.45 with a standard deviation of 0.37. Mean and standard deviation of complexity are 8.7 and 1.17, respectively. The creek

Geophysical Research Letters
networks demonstrate a clear bimodal distribution: network group I with higher complexity and efficiency and network group II with lower complexity and efficiency.
Two heat maps are created by summing the networks of each group to explore the probabilities of locations within the site being integrated into a creek network. Compared to group I, creek networks in group II are shorter and less complex (differences are mainly seen in areas A and B in Figure 2). An important location is identified and marked "C" in Figure 2 which acts as a threshold between groups I and II in area A. Specifically, simulated creeks that reach this point are then highly likely to continue expanding further inland, whereas those that do not, find no alternative route westward. In the area to the south of the wetland (area B), group II networks are restricted to topographic pathways close to the coast while group I networks extend across a greater area. From visual inspection, the bright areas in the heat maps show a high correlation with the existing drainage channels and low areas in the site. Figure 3 shows eight creek networks randomly selected from the 100 networks created by Case 2. Compared to those of Case 1, with the breach kept at the same location but without including preexisting channels in the initial site configuration, pathways of creek development are characterized by a higher level of meandering and branching. In this case, the mean of geometry efficiency is 2.74 with a standard deviation of 0.43, and the mean of complexity is 9.3 with a standard deviation of 0.83. Similar to the case mentioned above, creek networks in this case also show a clear bimodal distribution (Figure 4), with networks in group I more complex and efficient than those in group II. Remarkably, the two heat maps of the two network groups also identify location "C" as the threshold between network groups I and II in area A. The difference between the network groups in area B, that is, group I networks span across a larger area than group II networks, is also observed. Notably, the bright areas in the heat maps of this case also agree well with the existing drainage channels and low areas in the site. Figure 5 shows heat maps generated by summing the entire creek network ensemble of Cases 1 and 2. Although, as mentioned previously, pathways of creek development in the area covering the artificial short channels (zone 1) are very different in the two cases, relative to the effect in zone 1, it is clear that pathways of creek development in the rest of the floodplain are affected far less by the initial configuration of the site. The pattern in zone 2 is similar in both cases, that is, with or without the artificial channels adjacent to the breach. The corridor (feature 3) through which the creek network expands into the area to the south shows a high level of certainty. Bright edges of the wetland (feature 4) are also observed in both cases, outlining the topographic limit of the creek networks.

What Effect Does Vegetation Have?
Entire creek network ensembles of Cases 3 and 4, and the differences in probability of creek formation between Cases 3 and 1 and Cases 4 and 2 are shown in Figure 6. Pathways of creek development of Cases 3 and 4 are not dissimilar to those of Cases 1 and 2 ( Figure 5). However, probabilities of creek formation are higher in Cases 3 and 4. In both cases without vegetation (Cases 3 and 4), pronounced increases in probabilities are observed in a remote part of the network system (area A) beyond the above-mentioned location "C". Increased probabilities are also observed in area B. Further, locations with increased probabilities in area A follow the lines of the existing drainage channels. Geometry efficiencies and complexities of the 100 creek networks of Cases 3 and 4 are also given in Figures 6e and 6f. It is observed that when there is no vegetation, the creek networks are no longer separated into two groups. They show characteristics of the above-mentioned group I creeks-long and more complex.

Discussion and Conclusions
This research focuses on the development of a simple rules, reduced complexity creek evolution model that combines a computationally inexpensive hydrodynamic module, a full vegetation growth module, and a stochastic creek evolution module. The results illustrate the application of the model which emphasizes the use of Monte-Carlo simulation. The processes and phenomena the model aims to simulate, that is, the formation and evolution of creek networks, are characterized by high levels of uncertainty (Whitehouse et al., 2000;Dale et al., 2018). These uncertainties are accounted for through random headward expansion direction selection in each individual model run and represented as heat maps by way of Monte-Carlo simulation.
Despite the above-mentioned uncertainties, a certain relationship between tidal creek pathways and surface features is shown in our results, manifested by cases with different domain configurations producing heat maps with similar patterns. The bright areas in the heat maps correspond to the existing drainage channels and low areas in the site. These results indicate that although the processes the model aims to simulate are subject to random disturbances, there is a clear tendency towards specific tidal creek formation topologies with high probability determined by wetland topography.
Compared to wetland topography, the influence of vegetation is secondary in determining the heat map patterns. Vegetation, however, reduces the probabilities of tidal creek formation. In cases with vegetation, ensemble results identified a threshold point (location "C" mentioned above) in a remote part of the model domain which the simulated creek networks may or may not be able to reach and expand further inland. Further, in area B creek expansion is similarly impeded. Such weakened creek expansion due to vegetation has been observed in previous work (e.g., Wallace et al., 2005).
Although the distribution of vegetation responds directly to flow and inundation stresses which are affected by topography, distribution of vegetation is not directly constrained by topography, for example, salt marsh vegetation is present on higher ground. Spatial variation of vegetation stem density (see Figure S4) on the unchanneled platform then decreases rapidly over time due to lateral expansion. The initial spatial establishment of vegetation thus has a very small impact on the intermediate and final distribution of vegetation. Furthermore, resistance of the soil to erosion over the domain is increased in cases with vegetation. With this domain-wide increase in soil resistance, the threshold points which were previously (no vegetation) easily eroded are revealed (sometimes eroding and sometimes not eroding), leading to the bimodal creek formations.
Note that this vegetation implementation is not necessarily unrealistic because high ground, which might not be suitable for the simulated salt marsh vegetation to become established, could just as well be populated by other plants with lower salt tolerance that would affect the erodibility of the soil in a similar way. However, if the distribution of vegetation was further constrained by the topography, meaning more than one type of plant were grouped into patches in the domain, the impact of vegetation on the passing flow as well as on the erodibility of the soil (reflected by critical bed shear stress) should be differentiated. For example, fully submerged vegetation reduces the velocity of the passing flow and facilitates sedimentation, whereas semi-submerged vegetation also redirects flow and increases the probabilities of adjacent areas turning into creeks (Schwarz et al., 2014). Also, plants with shallow/sparse horizontal root branching may be less effective in reducing lateral/over edge erosion than those with deep/dense horizontal root branching (e.g., van Eerdt, 1985). By incorporating the plant model of Temmerman et al. (2007), the model could be applied to comprehensive vegetation-sensitive salt marsh ecosystem studies. For instance, realizing events that affect flow and inundation stresses, both in magnitude and duration (including storms), would enable consideration of climate change impacts.
An example of the practical significance of utilizing ensemble results from this stochastic simulation is seen in Figures 2 and 4. Although out of the scope of this paper, the grouped results, which would be unavailable through a deterministic simulation approach, could be used in numerous ways for decision making. For example, if a long-term objective was the formation of a complex, far-reaching tidal creek network (as in Figure 2. AI), a management scheme could be implemented to artificially assist creek formation. In this case the solution could be to clear a channel at point C or create channels in area B in Figures 2 and 4, as opposed to creating short channels near the breach. These could represent channel excavation plans for salt marsh restoration scenarios in which some parts of the sites are underlain by agricultural hard pans located very close to the surface, inhibiting channel development and rapid ecological recovery. Without these data, inducing the desired creek formations artificially would either have an undetermined (and likely lower) chance of success or come at greater financial cost.
The reduced complexity model developed herein can be operated by nonspecialist modelers (e.g., coastal managers) using routinely available data on wetland topography, vegetation distribution, and so on and thus can function as a management scoping tool. To further improve the model, future work that offers a promising avenue of investigation includes adding (1) more biological processes to the model, such as impacts of extracellular polymeric substances (EPS) on sediment erosion; and (2) sediment redistribution influence on the evolution of topography.