Modeling Oyster Reef Restoration: Larval Supply and Reef Geometry Jointly Determine Population Resilience and Performance

Restoration of native oyster (Crassostrea virginica) populations in Chesapeake Bay shows great promise after three decades of failed attempts. Population models used to inform oyster restoration had integrated reef habitat quality, demonstrating that reef height determines oyster population persistence and resilience. Larval recruitment drives population dynamics of marine species, yet its impact with reef height and sediment deposition upon reef restoration is unknown. To assess the influence of reef height, sediment deposition and larval supply, we adapted a single-stage population model to incorporate stage structure using a system of four differential equations modeling change in juvenile density (J), and changes in volume of adults (A), oyster shell reef (R), and sediment (S) on an oyster reef. The JARS model was parameterized with empirical data from field experiments. Larval supply included larvae from the natal population and from outside populations. The stage-structured model possessed multiple non-negative equilibria (i.e., alternative stable states). Different initial conditions (e.g., oyster shell reef height) resulted in different final states. The main novel findings were that the critical reef height for population persistence and resilience was jointly dependent on sediment input and larval supply. A critical minimum larval supply was necessary for a reef to persist, even when initial sediment deposition was zero. As larval supply increased, the initial reef height needed for reef persistence was lowered, and oyster reef resilience was enhanced. A restoration oyster reef with higher larval influx could recover from more severe disturbances than a reef with lower larval influx. To prevent local extinction and assure a positive population state, higher levels of larval supply were required at greater sediment concentrations to overcome the negative effects of sediment accumulation on the reef. In addition, reef persistence was negatively related to sediment deposited on a reef prior to larval settlement and recruitment, implying that restoration reefs should be constructed immediately before settlement and recruitment to minimize sediment accumulation on a reef before settlement. These findings are valuable in oyster reef restoration because they can guide reef construction relative to larval supply and sediment deposition on a reef to yield effective and cost-efficient restoration strategies.

Restoration of native oyster (Crassostrea virginica) populations in Chesapeake Bay shows great promise after three decades of failed attempts. Population models used to inform oyster restoration had integrated reef habitat quality, demonstrating that reef height determines oyster population persistence and resilience. Larval recruitment drives population dynamics of marine species, yet its impact with reef height and sediment deposition upon reef restoration is unknown. To assess the influence of reef height, sediment deposition and larval supply, we adapted a single-stage population model to incorporate stage structure using a system of four differential equations modeling change in juvenile density (J), and changes in volume of adults (A), oyster shell reef (R), and sediment (S) on an oyster reef. The JARS model was parameterized with empirical data from field experiments. Larval supply included larvae from the natal population and from outside populations. The stage-structured model possessed multiple non-negative equilibria (i.e., alternative stable states). Different initial conditions (e.g., oyster shell reef height) resulted in different final states. The main novel findings were that the critical reef height for population persistence and resilience was jointly dependent on sediment input and larval supply. A critical minimum larval supply was necessary for a reef to persist, even when initial sediment deposition was zero. As larval supply increased, the initial reef height needed for reef persistence was lowered, and oyster reef resilience was enhanced. A restoration oyster reef with higher larval influx could recover from more severe disturbances than a reef with lower larval influx. To prevent local extinction and assure a positive population state, higher levels of larval supply were required at greater sediment concentrations to overcome the negative effects of sediment accumulation on the reef. In addition, reef persistence was negatively related to sediment deposited on a reef prior to larval settlement and recruitment, implying that restoration reefs should be constructed immediately before settlement and recruitment to minimize sediment accumulation on a reef before settlement. These findings are valuable in oyster reef restoration because they can guide reef construction relative to larval supply and sediment deposition on a reef to yield effective and cost-efficient restoration strategies.

INTRODUCTION
Native oyster species are dominant ecosystem engineers, which provide a diverse suite of ecosystem services, including water filtration, food and habitat for various species, shoreline stabilization, coastal defense, and enhanced fisheries production (Cerco and Noel, 2007;Coen et al., 2007;Grabowski and Peterson, 2007;Grabowski et al., 2012). These species were once abundant worldwide, such as the eastern oyster Crassostrea virginica of the western Atlantic and Gulf of Mexico coasts (Winslow, 1881;Baylor, 1895). Unfortunately, native oysters have been nearly extirpated by overfishing, habitat degradation, disease and eutrophication (Rothschild et al., 1994;Jackson et al., 2001;Kirby, 2004;Lotze et al., 2006;Airoldi and Beck, 2007;Beck et al., 2011;Zu Ermgassen et al., 2012). Worldwide, 85% of native oyster reefs have been eradicated (Beck et al., 2011). Along the Pacific, Atlantic, and Gulf of Mexico coasts of North America, 88% of oyster biomass and 64% of reef area have been lost since the 1800s by the eastern oyster C. virginica and the Olympia oyster Ostrea lurida (Zu Ermgassen et al., 2012).
In Chesapeake Bay, abundance of the eastern oyster has been reduced to approximately 1% of its historical levels due to overfishing, oyster reef degradation, and disease (Rothschild et al., 1994;Wilberg et al., 2011;Schulte, 2017). Restoration efforts had failed through the start of the 21st century (Kennedy et al., 2011), but have proven successful since the mid-2000s due to the adoption of new restoration approaches in Chesapeake Bay (Schulte et al., 2009;Lipcius et al., 2015Lipcius et al., , 2019Colden et al., 2016;Theuerkauf and Lipcius, 2016;Lipcius and Burke, 2018) and in other Atlantic coast estuaries (Taylor and Bushek, 2008;Powers et al., 2009), and in the Gulf of Mexico (La Peyre et al., 2014). These efforts are distinguishable from past failed attempts by their attention to key elements of habitat quality, such as the architecture and height of restoration reefs (Taylor and Bushek, 2008;Powers et al., 2009;Schulte et al., 2009;La Peyre et al., 2014;Lipcius et al., 2015), and to insights gained from population models.
Several modeling approaches have been utilized to inform oyster population dynamics, restoration and fishery management strategies, each with specific goals and approaches. Mechanistic, physiologically-based models were among the first to be employed to address issues regarding effects of environmental conditions and biotic processes such as food availability (Hofmann et al., , 1994Powell et al., 1992;Dekshenieks et al., 1993) and later to the effects of disease and environmental factors (Powell et al., 1995(Powell et al., , 2011Dekshenieks et al., 1996;Klinck et al., 2002;Hofmann et al., 2004;Lavaud et al., 2017). Similar models were subsequently used to address management strategies and dynamics of shell budgets on harvested oyster reefs (Powell et al., 2006(Powell et al., , 2012(Powell et al., , 2018Powell and Klinck, 2007;Wilberg et al., 2013;Soniat et al., 2014). More recently, population models have been addressing issues related to oyster restoration, specifically the factors that promote long-term persistence and resilience of restoration reefs (Jordan-Cooley et al., 2011;Puckett and Eggleston, 2012;Pine et al., 2015;Housego and Rosman, 2016;Moore et al., 2016Moore et al., , 2018DePiper et al., 2017;Lipcius et al., 2019;Yurek et al., 2021). Some population models have integrated reef habitat quality (Jordan-Cooley et al., 2011;Wilberg et al., 2013;Housego and Rosman, 2016;DePiper et al., 2017;Lipcius et al., 2019;Yurek et al., 2021), and demonstrated that initial reef height is a key determinant of oyster population persistence and can lead to alternative states (i.e., reef degradation or persistence) for the eastern oyster in restoration (Jordan-Cooley et al., 2011;Housego and Rosman, 2016) and harvest (Wilberg et al., 2013;DePiper et al., 2017) reefs. Population dynamics of oysters and other invertebrates, such as corals, can differ as a function of age or stage (Moore et al., 2016(Moore et al., , 2018 and the magnitude of recruitment (Colden et al., 2017;Edmunds and Riegl, 2020). Hence, a need exists to incorporate stage or age structure and recruitment variation in models of eastern oyster population dynamics to determine optimal reef geometries in restoration. Consequently, we modified our previous model (Jordan-Cooley et al., 2011) by separating the population into juvenile and adult stages, which allowed us to assess the influence of variable recruitment and sediment concentration on critical reef height for population persistence of eastern oyster on restoration reefs. We hypothesized that the critical reef height for population persistence is not constant in relation to sediment being deposited on a reef, and is related to the magnitude of larval supply and recruitment impinging on an oyster reef.
To evaluate the hypotheses, we first describe how we integrated the juvenile and adult elements of the stage-structured model. Next, we define the parameters and functions, and show how they were derived. Finally, we conduct numerical and bifurcation analyses of the model for initial reef height under differing levels of larval supply and sediment concentration to produce a comprehensive evaluation of population dynamics relevant for oyster restoration.

STAGE-STRUCTURED DEMOGRAPHIC MODEL
Our original model consisted of three differential equations that describe the interactions between live oyster volume, dead oyster shell volume in the reef, and sediment accumulation (Jordan-Cooley et al., 2011). We now revise the model to incorporate stage structure by splitting the live oyster equation into two equations, one for adult volume and another for juvenile density. The revised JARS (Juveniles, Adults, Reef, Sediment) model defines the change in juvenile oyster density J, adult oyster volume A, dead oyster shell volume R, and sediment volume S on a 1-m 2 reef patch as a function of time t in years. The units for A, R and S may be interpreted as a volume (m 3 ) or as total reef height (m = m 3 · m −2 ) because we are dealing with a 1-m 2 reef patch. See Figure 1 for a schematic of the system and Table 1 for a summary of variables.

Juvenile Oyster Density
The rate of change in live juvenile oyster density J in number per m 2 is: FIGURE 1 | Representation of the main elements of the stage-structured JARS model; letters above each element match those in the model. Larval production P, larval supply, which includes larvae both from the natal population and from outside populations. C, water-column sediment concentration before being deposited on the reef. Photo by Stopher Slade; https://ncseagrant.ncsu.edu/coastwatch/previous-issues/2017-2/holiday-2017/making-the-most-of-oyster-reef-filtering/.
where the first term represents the product of larval production (= supply) P, a larval recruitment function L(A, R), and a switching function f (d J ). The second term (−mJ) defines a loss through maturation of juveniles into adults. Larval production P represents availability of larvae to settle onto a reef (Connolly and Roughgarden, 1999). In our model, larval production is not explicitly linked to adult abundance on the natal reef, such that larval production implicitly includes larvae supplied from outside reefs and the natal reef. The larval recruitment function L (i.e., the fraction of oyster larvae that can settle successfully, become juveniles, and survive to mature to the adult stage) is dependent upon the volume of live adult oysters A and the volume of dead oyster shell R, as: where ψ, ζ , and θ are fitting parameters (see Table 2 for list of parameters). In Figure 2, we display the form of the function L(A, R) with parameter values that best match field data. For a fixed R, larval recruitment with respect to A is a Ricker-type function (Schulte et al., 2009), such that recruitment peaks at some intermediate value of A and then decays as A → ∞ (Figure 2). Larval recruitment with respect to R, for fixed A, is a  sigmoid function (Colden et al., 2017), which increases slowly at low dead oyster shell volume, then rises sharply before saturating asymptotically as R → ∞ (Figure 2). The sigmoid curve of recruitment increasing toward an asymptote with only dead shell was based on the data in Colden et al. (2017) when only shell was present on restoration reefs in the first year of the experiment. The effects of adults were based on Schulte et al. (2009) and on the second-year results in Colden et al. (2017) when adults were present. Moreover, because the dependence of recruitment on adult oysters is a Ricker-type function, which is linear at smaller values of adults, and the dependence on dead shell is sigmoidal, then small to moderate numbers of adults are more favorable to recruitment than dead shell. This collective effect upon recruitment thus accounts for the preference of pediveligers to settle amongst conspecifics. When ζ > θ −1 θ , as is the case for our best fit parameters, L < 1 for any A, R ≥ 0 and L(A, R) can be interpreted as the fraction of available larvae that are successfully recruited and survive to adulthood. The parameter θ > 1 controls the behavior of decay with A; larger θ produces faster decay. The parameter ζ > 0 determines the height of the peak, while ψ > 0 defines the location of the peak. The collective effect of both functions is shown in Figure 2. The maturation term (−mJ) in Equation (1) contains a maturation rate m. Because the input term for the juvenile stage is larvae that are successfully recruited and will survive to adulthood, all of the J stage transitions to the adult stage.

Switching Functions
The proportion of oysters unaffected by sediment, which determines an effective population growth rate of juvenile or adult oysters, is given by a switching function: where d i = a i A + R − S and i = J, A. Here d J and d A represent the volume of live oysters and dead oyster shell not affected by Frontiers in Marine Science | www.frontiersin.org  Table 2. sediment on the reef, and a A > a J . The proportion of oysters not affected by sediment is an increasing function of d A or d J with a sigmoid shape bounded by 0 and 1 (Figure 3). When deposited sediment covers both the dead oyster shell (R) and a fraction a i of the adult oyster volume (a i A), d i is near 0. When d i > 0, live oysters are not significantly covered by sediment. The value of f (d i ) approaches 1 and the oysters can grow. In contrast, if more than a fraction a i of the adult oysters are covered by sediment, then d i < 0, f (d i ) approaches 0 and the oysters cannot grow. The same switching function was also used in Jordan-Cooley et al. (2011), but here the switching function for juveniles and adults differs in the parameters a J and a A , which set the burial threshold for detrimental effects of sediment on juveniles and adults, respectively. Specifically, the threshold for juveniles is much lower than that for adults ( Table 2).

Live Adult Oyster Volume
The change in live adult oyster volume A is described by: The first term (αmJ) converts the number of matured juveniles to volume of adult oysters. The second term is the population growth rate of adult oysters, which like in the previous model (Jordan-Cooley et al., 2011) assumes live oysters follow logistic growth with an intrinsic growth rate φ and carrying capacity K. We call the growth rate parameter φ instead of r here because it represents only the positive change in biomass, and does not include the per capita mortality rate µ due to predation and disease incorporated in r (r = φ − µ). Population growth and mortality are scaled by a factor f (d A ), which is defined in Equation (3), since only oysters above the sediment can survive and reproduce. The last term represents the decrease in live adult oyster volume due to suffocation by sediments, where ǫ is the mortality rate of oysters covered by sediment. This term has a stronger effect as f (d A ) becomes smaller; i.e., when more oysters are covered in sediment.

Dead Oyster Shell Volume
The change of dead oyster shell volume R is modeled as:  Table 2.
where the first two terms in the equation are directly from mortality of live oysters. The third term is the loss of dead oyster shell volume due to degradation, which is proportional to the volume of dead oysters at the rate γ . For simplicity, we do not include juveniles in the dead shell volume because they are small enough that they do not contribute significantly to shell volume .

Sediment Accumulation
The change in sediment volume S is: where the first term represents the rate of sediment erosion, which is proportional to the volume of sediment on the reef at a rate β. The second term is the rate of sediment deposition. Without filtration by live oysters, the sediment deposition rate is Cg, where C is a maximum possible deposition rate and g is a modification function that decreases as reef height (A + R) increases. The rate of sediment deposition is maximized when there is no reef and no oysters, in which case the deposition rate is just C. Typical values of sediment deposition range from 0.00 to 0.02 m y −1 (Cronin et al., 2003), which is why we chose a maximum rate for C of 0.02 m y −1 ( Table 2). The parameter C is the sediment deposition rate from the water column (Figure 1), whereas the amount of sediment actually deposited and remaining on the reef is the variable S, which is a non-linear function of C. The effect of sediment on larval settlement and recruitment is implemented via the switching function f (d i ) in the input term of the juvenile equation. When the sediment height is small relative to the reef height, it is assumed that sediment has minimal effect on the larvae. However, f (d i ) causes the fraction of available larvae successfully recruited to decline as sediment increases, reflecting the realistic situation in nature.
A filtration rate function F scales down the deposition rate due to filtration of sediment by adult oysters A. The filtration function F and modification function g are defined by The higher the filtration rate, the smaller the deposition rate. The filtration rate F scales linearly with Cg when Cg is small, reaches a peak filtration rate F 0 at some optimal sediment concentration y 0 , and beyond this threshold, it decreases as oyster gills become increasingly clogged (Jordan-Cooley et al., 2011). The premise is to have a function bounded between 0 and 1 that peaks at some intermediate value and drops off exponentially. Graphs of g and F(Cg) are displayed in Jordan-Cooley et al. (2011).
In the Appendix, we show that a solution of Equations (1), (4), (5), and (6) with non-negative initial values is always non-negative and bounded; hence the model is mathematically well-posed.

PARAMETER ESTIMATES
In the revised JARS model (Equations 1, 4, 5 and 6), some parameters come directly from field observations or the literature ( Table 2), whereas others are approximated by fitting the functions with data from field restoration reefs, as described below. Although the model parameters were estimated from a limited number of sites in Chesapeake Bay, the selected sites are the best-studied, long-term restoration sites in the bay (Schulte et al., 2009;Lipcius et al., 2015Lipcius et al., , 2019Colden et al., 2016;Theuerkauf and Lipcius, 2016).

Maturation Rate
We assume the average maturation time to adulthood is 1 year, such that the parameter m = 1 yr −1 . The maturation rate does not affect the equilibrium values of adults, reef or sediment because the input term in the juvenile equation is the rate at which juveniles are successfully recruited. All juveniles will eventually move into the adult stage regardless of the speed at which they do so, because the juvenile equation does not include mortality. Changing the maturation rate would only change the number of juveniles in our equilibrium results. However, it is conceivable that changing the maturation rate could affect transient dynamics. Using an extreme maturation rate of m = 3 yr −1 , we reran our calculation of the needed height of dead shell for successful restoration, and the results were virtually identical to those presented here.

Conversion of Juvenile Density to Oyster Volume
Since the units for juvenile oysters J and adult oysters A are numbers · m −2 and m 3 · m −2 respectively, we needed a conversion factor α in Equation (4) to convert juvenile density into adult volume. To do so, we used shell height and volume data for 28 oysters from the field study of Lipcius et al. (2015). The data were log 10 -transformed and analyzed with linear least-squares regression (Figure 4), yielding the following equation: where V is shell volume in mm 3 and SH is Shell Height in mm. From this relationship, we calculated the volume of one juvenile oyster using its shell height. We assumed that mean juvenile SH = 35 mm (Schulte et al., 2009). Using Equation (8), the volume of an average juvenile oyster is 3.256 · 35 2.184 · 10 −9 = 7.6724 · 10 −6 m 3 . When converting, we also needed to multiply the calculated juvenile volume by a packing parameter p (Ducharme-Barth et al., unpublished manuscript) 1 , because the oysters are not all packed together and there is interstitial space between individuals. Using the packing parameter 4.6 (Ducharme-Barth et al., unpublished manuscript) 1 , we scaled the volume with α = 7.6724 · 10 −6 · 4.6 = 3.529 · 10 −5 m 3 per juvenile.

Larval Recruitment Function
To estimate fixed parameter values for ψ, θ , and ζ in the larval recruitment function L(A, R) and site-dependent values for larval availability P (Equation 2), we used a Matlab function to fit data collected at different restoration sites in the Great Wicomico River (GWR) by Schulte et al. (2009), and in the Great Wicomico River, Linkhorn Bay (LB) and Broad Bay (BB) of the Lynnhaven River system by Colden et al. (2017). Our approach was to fit the predicted number of recruits PL(A, R) in each location and year to the observed number. Because successful recruitment occurred, we assumed sediment did not inhibit recruitment and thus ignored the switching function f (d) in the recruitment term of Equation (1). Available data for adults were number · m −2 , which we converted to m 3 · m −2 as described above and assuming a mean adult shell height of 100 mm Lipcius and Burke, 2018). One of the sites (GWR3) did not have data available for dead shell volume, so we assumed the shell volume R was 25% of the adult volume (Schulte et al., 2009). We obtained similar results when that data set was excluded.
Now having values for A and R in m 3 · m −2 and for J in number· m −2 , we imported the data into Matlab and sorted the data with respect to location and whether adults were present. Then, we defined a cost function : where the purpose of is to capture squared differences between model and data. To estimate recruitment without oysters, we used data from Colden et al. (2017) for the first year of the study when only shell was present at the start of the experiment.
To estimate recruitment with adults present, we used data from Schulte et al. (2009) and Colden et al. (2017) after the first year when adults were present on the reefs. In the equation, i values from 1 to 7 are the sites 1 (GWR3 with shell and adults), 2 (GWR1 and GWR2 with shell and adults), 3 (LB with shell and adults), 4 (BB with shell and adults), 5 (GWR1 and GWR2 with shell, without adults), 6 (LB with shell, without adults), and 7 (BB with shell, without adults); j is the index of data points in each data set and k i is the number of data points in the ith data set. Data

GWR3 GWR1 and GWR2 LB BB
With adults P(1) = 570 P(2) = 112 P(3) = 326 P(4) = 238 Without adults −− P(5) = 1, 086 P(6) = 513 P(7) = 1, 884 from all sites were combined into a single data set for fitting. Because GWR3 produced more data than the other locations, we also conducted the analysis without GWR3 data for comparison. Estimated parameter values for L(A, R) were similar when the GWR3 data were excluded, as were critical reef heights. The data sets satisfy k 1 = 84, k 2 = k 5 = 12, and k 3 = k 4 = k 6 = k 7 = 6. For each data index ij, P(i) · L(A ij , R ij ) is predicted juvenile density after a year of recruitment (calculated using the model parameters and reef structure at site i), and J ij is the actual data point representing juvenile density. The exponential terms were added in Equation (9) to ensure that all parameters in L(A, R) estimated in the optimization process were positive. The data sets with values of (A ij , R ij , J ij ) for 1 ≤ i ≤ 4 and (0, R ij , J ij ) for 5 ≤ i ≤ 7 were input into and a Matlab function fminsearch was used to minimize the function to find the best fitting parameters (ψ, θ , ζ ) and P(i) with 1 ≤ i ≤ 7. Parameter estimates of the fitted functions differed significantly from zero. The fitted function for L(A, R) produced estimates of ψ = 0.3135, θ = 2.4300, and ζ = 0.8068. Estimates for P are in Table 3. Fitting results are shown in Figure 2 where L(A, R) is plotted as a function of both A and R. To illustrate how this compares with the data, the observed and fitted yearly recruitment at each site are shown in Figure 5 (plotted vs. the combination A + R for convenience of visualization).

Carrying Capacity
We estimated carrying capacity K by assuming that mean SH = 100 mm for adult oysters (Schulte et al., 2009). Using the conversion function (Equation 8), the volume of an average adult oyster = 3.256 · 100 2.184 · 10 −9 = 7.5977 · 10 −5 m 3 . From field data, the maximum number of adult oysters in 1 m 2 averages 200 (Schulte et al., 2009). Including the packing parameter, K = 7.5977 · 10 −5 · 4.6 · 200 = 0.06990 m 3 · m −2 , which we rounded up to 0.1 as a maximum volume for K. We note that because K appears in the term of Equation (4) for positive change in biomass, it is not a carrying capacity in the traditional sense that populations will approach this value after a long time. Rather, it represents the population above which biomass cannot be generated.

Dead Oyster Shell Degradation Rate
For the dead oyster shell volume equation (Equation 5), the oyster shell degradation rate γ is estimated with the assumption that the average half life of oyster shell is 6 y (Powell et al., 2006). In the absence of live oysters, Equation (5) becomes a linear decay equation dR dt = −γ R. Solving it with a half life of t = 6 y yields R(t) = R(0)e −γ t = R(0)/2, such that γ = ln(2)/6 = 0.1155 y −1 .

Consistency of Single-Stage and Stage-Structured Models
Similar to the original model (Jordan-Cooley et al., 2011), the stage-structured JARS model possessed multiple equilibria. The extinction equilibrium E 0 = (0, 0, 0, C/β) always existed, and for most reasonable parameter values, it was locally asymptotically stable (see Appendix). This implies that for initial values close to the extinction state, the population would always go extinct. Moreover, persistence of the population occurred under similar values of C and initial conditions as in Jordan-Cooley et al. (2011). Below we document the regimes for parameters P (larval production) and C (water-column sediment deposition rate) that supported conditional persistence of the population, thus inducing alternative stable states. As with C, the parameter P was chosen due to its variability among oyster populations. When alternative stable states existed, we determined the initial conditions that led to persistence or extinction. While our results may also be applicable to a wider range of parameter values, we used the parameter values in Table 2 because they are based on previous studies, field data, and our estimates from fitting field data for restoration reefs.

Bifurcation Analysis
The effects of parameters P and C on the oyster population were examined through bifurcation analysis to determine fourdimensional (J, A, R, S) equilibrium population sizes as a function of the bifurcation parameters P and C. Bifurcation analysis employs bifurcation theory in non-linear dynamical systems (Kuznetsov, 1998) to produce a comprehensive evaluation of population dynamics. In population models, bifurcation analysis can detect ecological steady states as a function of model parameters and bifurcation points where dynamics change drastically (de Roos, 2021). Bifurcation diagrams were generated by Matlab package Matcont (Dhooge et al., 2003).
A bifurcation diagram with respect to larval production P is shown in Figure 6. In this and subsequent bifurcation diagrams, equilibrium population sizes (solid stable and dashed unstable curves) are depicted as a function of specific bifurcation parameters. Each bifurcation diagram also displays a saddlenode bifurcation point (LP) that represents the boundary between alternative stable states and population extinction (i.e., the tipping point). For instance, in Figure 6, P is the bifurcation parameter, such that for each value of P, equilibrium population sizes are indicated by the curves. The equilibria are four-dimensional (J, A, R, S), so we display each of the variables separately for clarity. Resilience is defined as the ability of a population to recover from a disturbance that reduces population abundance significantly. In Figure 6 and other bifurcation diagrams, resilience is proportional to the distance between the stable and unstable positive equilibria (Scheffer and Carpenter, 2003).

Bifurcation Structure as a Function of P
The parameter P in Equation (1) represents the rate of larval production, and scales the larval recruitment function L(A, R). For P there was a threshold value P critical = 220.6 such that when P < P critical , there was no equilibrium with a positive population size (Figure 6); i.e., insufficient numbers of larvae were available to settle either from the natal population or from outside populations. All positive initial conditions approached an extinct state with zero population and high sediment level (not shown in the graphs). When P > P critical , the system had two positive equilibria (Figure 6). The upper equilibrium (solid blue curve) represented stable high population sizes to which sufficiently large initial population sizes were attracted. The lower equilibrium (dashed orange curve) was an unstable smaller population representing a threshold below which populations went to extinction; above it populations moved toward the stable positive equilibrium (blue curve). Also, resilience increased proportional to larval production P (Figure 6), as reflected in the greater distances between stable solid blue and unstable dashed orange curves at high P. P(i) values obtained through fitting field data in Table 3 were all above this threshold value P critical = 220.6, except that P(2) was below. For P = 250, which will be used later as a benchmark point, there were two positive equilibrium states: (J, A, R, S) = (93.41, 0.05, 0.16, 0.02), which is locally asymptotically stable, and (J, A, R, S) = (73.8, 0.04, 0.15, 0.05), which is unstable. Similar bifurcation structure resulted when we used different parameter values in the analysis.

Bifurcation Structure as a Function of C and P
In our original single-stage model (Jordan-Cooley et al., 2011), the bifurcation structure of positive steady states of the model was explored in terms of maximum sediment deposition rate C. The bifurcation structure with parameter C was preserved in the JARS model (Figure 7). There was a threshold value C critical , which depended on P, such that when C < C critical , the system had two positive equilibria and resilience of J, A, and R increased as C decreased (Figure 7). When C > C critical , there was no positive equilibrium; the extinction state attracted all initial values (Figure 7). The effect of the larval production rate P on bifurcation structure of C was to shift the saddlenode bifurcation point and positive equilibrium curves; larger values of P allowed the population to persist at higher values of C and increased resilience. For example, at P = 250, C critical = 0.0225 (Figures 7A-D), whereas at P = 800, C critical = 0.070 (Figures 7E-H). This level of larval production P = 800 is within the observed range of values in Chesapeake Bay (Schulte et al., 2009;Lipcius et al., 2015;Colden et al., 2017). The structure of the reef (values of J, A, and R) was insensitive to the sediment deposition rate C as long as C was below the bifurcation point (Figure 7). In contrast, the reef structure depended strongly on the larval production rate P.
Investigating the joint effects of P and C upon bifurcation structure also allowed us to obtain an estimate for the instantaneous growth rate of adult oysters φ. The threshold production rate of oyster larvae P for positive population states to exist was a function of maximum sediment deposition rate C (blue curve in Figure 9). To prevent local extinction and assure a positive equilibrium state, higher levels of P were required at higher concentrations of C to overcome the negative effects of sediment accumulation on the reef. Even at low maximum sediment deposition rates C, the larval production rate P had to be positive for populations to persist because the production of oyster larvae was the only source of new recruits in the model. Hence, we chose the instantaneous growth rate of adult oysters φ in Table 2 to reflect this requirement.
Although slightly non-linear, the location of the saddle-node bifurcation point (blue curve in Figure 9) was well-approximated by P = 11, 205 C. For this set of environmental conditions, this provided a criterion for oyster restoration: If P > 11, 205 C at a location, it would be possible to build a persisting oyster reef.

Threshold Reef Height
For application of the JARS model to oyster restoration, we determined critical values of R 0 , initial reef height of dead oyster shell, to construct a reef for successful restoration. Threshold R 0 values were computed numerically, assuming the parameter values in Table 2. When S 0 = 0, there is initially no sediment on the restoration reef (Figure 8). As an illustration of the numerical analysis used to determine threshold reef height R 0 , in Figure 8 we chose larval availability P to be 250, which is a scenario with sustainable dynamics (above P critical = 220.6). When R 0 = 0.26 m the population persisted (Figure 8A), whereas when R 0 = 0.25 m the population declined to extinction ( Figure 8B). Thus, the critical value should be between 0.25 and 0.26 m. The calculated critical value for R 0 was approximately R 0c = 0.253. In general and for other values of P, when R 0 ≥ R 0c (Figure 8A), the oyster population and reef persisted ( Figure 8A). After some time, it approached the stable positive state. Otherwise, sediment S accumulated rapidly and the system moved toward an extinct state ( Figure 8B).
We also evaluated the joint influence of maximum sediment deposition rate C and production of oyster larvae P on R 0c (Figure 9) because these two parameters can differ substantially between locations (Cronin et al., 2003;Lipcius et al., 2015;Colden et al., 2017). The main result from this simulation was that at a given level of sediment deposition rate C, there existed a wide range of R 0c depending on larval availability P. The higher the P, the lower the initial reef height needed for any level of C (Figure 9). For instance, at a maximum sediment deposition rate C = 0.02 and P = 250, the initial dead shell reef height R 0 needed for persistence was about 0.253 m, whereas with P = 1,000, R 0c  Table 2. Blue solid branch: stable; orange dashed branch: unstable.  Table 2. FIGURE 9 | Critical reef height R 0c as a function of C and P in a 2-d color gradient representation. Blue curve is a superimposed bifurcation curve showing the saddle-node bifurcation point with respect to the joint effects of parameters C and P. Non-extinct outcomes exist only above the curve (including in the region of very small sediment deposition rates where R 0c was not computed). Other parameters are as in Table 2. The vertical bar legend to the right of the graph represents a gradient of R 0 in m.  Table 2. = 0.128 m, which is half that needed for P = 250. However, at each level of P there existed a bifurcation point C critical such that at values of C above C critical the population moved toward extinction irrespective of R 0 (Figure 9). The results of Figure 9 were robust to different values for maturation rate m and different parameters in the larval recruitment function L(A, R) obtained from fitting the recruitment data with and without GWR3 data.
For oyster restoration, we also focused on defining critical values of R 0 , initial reef height of dead oyster shell, given different values of S 0 , the initial height of sediment on a reef before larval settlement and recruitment occur (Figure 10). Under this scenario, restoration reefs are constructed with shell but without live oysters. Differing values of S 0 are applicable when restoration reefs are constructed at different time periods prior to settlement and recruitment, and sediment deposition C results in accumulation of sediment S on the reef. In Figure 10A, we used the lowest possible value for larval influx P = 221 for which reefs would persist. In this case the initial dead shell reef height R 0 necessary for persistence increased rapidly from under 0.3 m at S 0 = 0 to nearly 0.4 m at S 0 = 0.03 m. At higher levels of larval influx, such as P = 1, 000 (Figure 10B), the critical reef height dropped significantly; e.g., at S 0 = 0 critical reef height was less than 0.13 m for P = 1, 000, which was less than half the 0.3 m needed for P = 221 (Figure 10). For all levels of P, critical reef height increased with initial sediment. Especially when larval supply was poor, much higher initial reef heights R 0 were required for persistence to counteract the effect of sediment buildup on the reef before settlement occurred, as the initial sediment volume S 0 was amplified (Figure 10).

DISCUSSION
A stage-structured model composed of juvenile density J, and volumes of adults A (i.e., surrogate for biomass), dead oyster shell reef R, and sediment S on an oyster reef (JARS model) was parameterized with empirical data from field experiments. This enabled us to determine the collective influence of larval supply P and sediment deposition C on oyster reef resilience, and the critical initial dead-shell reef height R 0 for reef persistence. Larval supply was assumed to include the collective supply from the natal population and from larval subsidies from outside populations. Resilience was defined as by Scheffer and Carpenter (2003), specifically the ability of the oyster reef population to recover from disturbances that significantly reduce population abundance.
Similar to the single-stage model (Jordan-Cooley et al., 2011), the stage-structured JARS model possessed multiple equilibria (i.e., alternative stable states), such that different initial conditions of dead-shell reef height produced different final states. The main novel findings from the JARS model were that the critical initial reef height R 0 for population persistence was dependent not only on sediment input C, as in Jordan-Cooley et al. (2011), but also crucially on larval supply P and on S 0 , the amount of sediment built up on a reef before settlement and recruitment of oysters. As larval influx P to a reef increased, the initial reef height R 0 needed for reef persistence was lower and oyster reef resilience was enhanced, such that an oyster reef with higher larval influx could recover from more severe disturbances than a reef with lower larval influx. In addition, a threshold larval supply P was identifiable, below which a shell reef could not survive. These results are consistent with field experiments using high-and low-relief reefs; the critical reef heights determined here were in the same range (0.15-0.30 m) as those in field experiments where reef height, larval recruitment and sediment concentration varied (Schulte et al., 2009;Lipcius et al., 2015;Colden et al., 2017). The effects of the magnitude of larval supply P on the impact of other parameters was to lower the magnitude of positive parameters, such as the instantaneous growth rate of adult oysters φ, required for population persistence. For example, larger values of P allowed the population to persist at lower values of φ and at higher values of negative parameters such as sediment concentration C, by increasing population resilience to disturbance. These findings are valuable for our understanding not only of the ecology of oyster populations, but especially for oyster reef restoration because they can guide reef construction and placement to be effective and cost-efficient. When R 0 was greater than the threshold value of R, the oyster population and reef persisted. Otherwise, sediment S accumulated rapidly and the reef degraded to extinction. For larger values of P, the critical reef height decreased at most sediment concentrations. For instance, the critical reef height with no sediment at P = 300 was 20% lower than that with P = 221, which could reduce restoration costs significantly. This would be the situation when restoration reefs are constructed only out of dead shell R with no adults A 0 = 0 and immediately before the settlement and recruitment season such that S 0 = 0. In contrast, if reefs are constructed well before the settlement and recruitment period, and sediment S is allowed to build up on a reef before settlement and recruitment, then the relationship between R 0 and S 0 becomes critical in defining the additional initial dead shell reef height R 0 necessary for reef persistence. In our analyses, the critical R 0 was positively related to S 0 such that much higher initial reef heights R 0 were needed for reef persistence when the initial sediment volume S 0 on a reef was substantially greater than 0.
The characteristics of the interactions observed in the JARS model are highly relevant for marine species exhibiting metapopulation source-sink dynamics, which are common in marine and estuarine ecosystems (Lipcius and Ralph, 2011). For instance, the Lynnhaven River in lower Chesapeake Bay harbors an oyster metapopulation consisting of populations in interconnected bays and tributaries . Connectivity among the populations, as determined by biophysical modeling with historical oyster reef information, was complex and indicative of metapopulation structure . The metapopulation was apparently composed of four types of reef subpopulations with different implications for larval exchange and consequences for restoration strategies. Source subpopulation reefs were capable of self-replenishment and were interconnected with many of the reefs in the subestuary via larval exchange, indicating that these reefs were not only self-replenishing but also supplied larval subsidies to other reefs. These reefs are optimal for habitat restoration R and broodstock A enhancement, especially where no oyster population exists and the metapopulation must be jump-started with new source subpopulations. In contrast, sink subpopulation reefs did not self-replenish and larvae from sink reefs were advected out of the system, which makes these reefs unsuitable for broodstock A enhancement. However, sink reefs are suitable for habitat R restoration because these reefs would receive larval subsidies from source reefs once the source subpopulation reefs become established. Other self-replenishing subpopulation reefs supplied adequate larvae back to the natal reef, but did not supply larvae to other reefs in the system. Although such reefs are acceptable for habitat restoration and broodstock enhancement, they are, however, suboptimal in that their value to the metapopulation is limited. Finally, putative source subpopulation reefs alternated between self-replenishment and provision of larvae to other reefs, depending on environmental conditions. These reefs are thus acceptable for restoration and broodstock enhancement, but only after source subpopulation reefs have been restored.
The implications of our model for restoration are consistent with field experiments using high-and low-relief reefs, and the critical reef height determined here was very close to the field results (Schulte et al., 2009;Lipcius et al., 2015;Colden et al., 2017). However, initial conditions could vary from one site to another, and reliable estimates of parameters such as larval supply and sediment deposition rate are needed to determine the corresponding critical reef heights. In addition, habitat quality must be incorporated through habitat suitability index models (Theuerkauf and Lipcius, 2016) to define the optimal restoration strategy. By taking into account oyster juveniles as well as larval recruitment, our JARS model incorporates the age structure of oyster populations and thus provides a better conceptual framework for real-world oyster restoration.
Future refinements of the model should include timevarying parameters. For instance, larval supply P is actually time varying due to the effects of salinity and hydrodynamic changes resulting from freshwater discharge and rainfall, so fluctuations in P should be investigated. Also, we didn't model the effect of non-zero initial adult oysters A on critical reef R height, which may be significant due to the filtering ability and reproductive output of adult oysters (Wilson, 2019). Furthermore, larval supply P could be separated into natal and outside components to reflect the likely different contributions of outside and natal populations to reef persistence. Such modifications could be effectively modeled using the JARS model system.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
RL provided field data used to parameterize the model and produced the function of oyster biomass vs. size. YZ, JZ, LS, and JS analyzed the model. RL, YZ, LS, and JS wrote the manuscript. All authors developed the model.