Onset of runaway fragmentation of salt marshes

Salt marshes are valuable but vulnerable coastal ecosystems that adapt to relative sea level rise (RSLR) by accumulating organic matter and inorganic sediment. The natural limit of these processes defines a threshold rate of RSLR beyond which marshes drown, resulting in ponding and conversion to open waters. We develop a simplified formulation for sediment transport across marshes to show that pond formation leads to runaway marsh fragmentation, a process characterized by a self-similar hierarchy of pond sizes with power-law distributions. We find the threshold for marsh fragmentation scales primarily with tidal range and that sediment supply is only relevant where tides are sufficient to transport sediment to the marsh interior. Thus the RSLR threshold is controlled by organic accretion in microtidal marshes regardless of the suspended sediment concentration at marsh edge. This explains the observed fragmentation of microtidal marshes and suggests a tipping point for widespread marsh loss.


Introduction
There is a growing consensus that marsh vulnerability to rela-2 tive sea level rise (RSLR) is tied to inorganic sediment availabil-3 ity 1-4 , where deposition of inorganic sediment increases with 4 flooding duration, and potentially offsets sea level rise. Indeed, 5 inorganic deposition rates have accelerated over the last cen- 6 tury concomitant with sea level rise 5,6 and historic marsh loss 7 has been observed (and projected 7,8 ) mostly in sediment-poor 8 systems 9, 10 and microtidal marshes 11 . Modeled threshold 9 rates of RSLR for marsh drowning, using simplified point (0-D) 10 models, increase by 2 orders of magnitude as a function of 11 suspended sediment concentration and tidal range 12,13 . How- 12 ever, a contrasting body of work emphasizes the importance 13 of organic matter accumulation in building marsh soils in the 14 face of sea level rise, especially in the sediment deficient estu-15 aries most vulnerable to sea level rise 1,11,[14][15][16][17] . Total marsh 16 accretion rates are more strongly correlated with the organic 17 fraction of marsh soil than the inorganic fraction 14 ; organic 18 matter contributes 4 times more soil volume than an equiva-19 lent mass of inorganic sediment 16 ; and organic matter is the 20 dominant contribution to marsh accretion by volume in many 21 Atlantic and Gulf Coast marshes [14][15][16] . 22 Competing ideas about the relative importance of organic 23 and inorganic accretion likely reflect strong spatial gradients 24 within marshes [18][19][20] . Inorganic accretion increases with sus- submergence that are sometimes difficult to explain. For ex-36 ample, marshes along the Blackwater River (MD, USA) are 37 rapidly submerging despite having a higher suspended sed-38 iment concentrations measured in channels, than in nearby 39 stable marshland 32, 33 . Elsewhere, marshes are submerging 40 despite measured accretion rates that are similar to or ex-41 ceed RSLR 2,33 , which suggests measurements take place 42 mostly along marsh edges, where maximum accretion rates 43 are generally observed 21, 23, 34, 35 . 44 The complexity of organic and inorganic accretion in a 45 marsh platform leads to the simple question: where in a marsh 46 should organic and inorganic contributions to marsh accre-47 tion be characterized to best evaluate marsh vulnerability to 48 RSLR? Measurements from high elevation portions of a marsh 49 potentially underestimate future marsh accretion because inor-50 ganic accretion rates may accelerate with increased flooding 51 duration 2 . However, if low elevation marshes are also closest 52 to channels, then accretion rates from low elevation portions 53 of the marsh would overestimate accretion to the marsh as a 54 whole, and lead to an underestimation of marsh vulnerability 55 to RSLR. 56 Another issue with the interpretation of measured accretion 57 rates is that they tend to converge towards the local rate of 58 RSLR, as the marsh platform approaches an equilibrium eleva-59 tion 36 , which complicates the estimation of maximum accretion 60 rates unless marshes are already drowning 2,37 . Thus, there 61 is a need for better numerical models that resolve the spatial 62 complexity of marsh sediment dynamics 4,13,19,27,28,[38][39][40] . 63 A few existing process-based models (e.g. 19, 28 ) capture the 64 observed drowning of interior marshes and their conversion to 65 ponds [41][42][43] . They suggest marsh drowning, and subsequent 66 The current understanding of the onset of marsh loss is 92 that it takes place whenever marsh depth relative to mean 93 high water is higher than a critical value D c above which  Table S1 for details) 42, 43, 46-48 . 101 Assuming the existence of a critical depth for marsh recov- In what follows we present and validate an explicit expres-122 sion for the inorganic accretion rate across the marsh platform 123 and use it to obtain the critical inorganic accretion rate for 124 marsh drowning. We then introduce the drowning threshold, 125 characterize the runaway marsh fragmentation regime and 126 discuss the effect of external parameters on marsh drowning. 127

128
Exponential decay of sediment concentration 129 As inorganic sediments in the water column settle on the marsh 130 surface, where erosion is assumed to be negligible 27 , the av-131 eraged sediment concentration C decays with the distance 132 x from the channel or tidal flat (Fig. 2). Sediment concen-133 tration thus reaches its lowest value at the location furthest 134 away-a distance L-from marsh edges ( Fig. 2A), defined in 135 the model as the watershed divide. This decay in sediment 136 concentration is well approximated by an exponential function, 137 C(x) = C(0)e x/L c (as proposed by 25 and observed by 23 ), with 138 decay length L c (see Experimental Procedures). Therefore, 139 the inorganic accretion rate for a non-flat marsh platform can 140 be approximated as where the average sediment concentration C(0) at the channel 142 bank or marsh edge is proportional to the average concentra-143 tion C 0 at the channel or mud flat during flood (see Fig. S3 for 144 the proportionality factor). 145 The decay length L c of the average suspended sediment 146 concentration scales as the ratio of the tidal discharge per unit 147 width and the effective sediment settling velocity w f , in agree-148 ment with the scaling of the deposition length in unidirectional 149 turbulent suspensions 50 (Experimental Procedures). We find 150 tidal discharge per unit width scales as Ld z/T , where d z is 151 the tidal range, T is the tidal period and L is the characteristic 152 length of the local drainage basin. Thus, the decay length has 153 the form with fitting parameter b ⇡ 1.5, in agreement with both numeri-155 cal simulations and analytical approximations (Experimental 156 Procedures and Fig. 2B).
We find the exponential approximation accurately describes 158 the sediment concentration profile except in the region around 159 the watershed divide, where tidal flow stops and the simulated 160 average sediment concentration, and thus accretion rates, 161 converge to zero (Fig. 2). In reality, complex tidal flows may 162 lead to residual accretion rates in the marsh interior (e.g. 22 ), in 163 which case the exponential approximation provides an upper 164 limit to evaluate the resiliency of drowning marshes. In what 165 follows we use the watershed divide as a formal definition of 166 the marsh interior. 167 The exponential decay correctly predicts the spatial gradient 168 in the average sediment concentration and inorganic accretion 169 rates for a wide variety of salt marshes (Fig. 3), including low-170 elevation micro-tidal marshes in the Virginia eastern shore 171 (Phillips Creek) 34 and Georgia 35 , and meso-and macro-tidal 172 marshes in Plum Island, MA 51 , Norfolk, UK 21 and in the Bay of 173 Fundy, CA 52 (see Experimental Procedures for further details 174 on the analysis and interpretation of inorganic accretion data). 175 The scaling of L c with the tidal range d z (Eq. 2) means that 176 suspended sediments deposit closer to channels (or tidal flats) 177 at lower tidal ranges, whereas they are more homogeneously 178 distributed at higher tidal ranges. This is consistent with the 179 2/15 trend observed in field measurements (Fig. 3)  length we get for the critical inorganic accretion rate: where`c = L c /L = b /w + f is the rescaled decay length, which 198 only depends on the rescaled effective falling velocity with t c ⌘ t(D c ). Thus, the critical inorganic accretion rate fragmentation regime strongly depends on whether perma-232 nent ponds are isolated or connected to the channel network 233 (Fig. 5A). In the first case, tidal basins and watershed divides 234 remain unchanged and the system evolves towards a new 235 equilibrium state (Fig. 5A, left). The portion of the marsh 236 closer to the edge adapts to RSLR and reaches a non-uniform 237 equilibrium marsh elevation in response to spatial gradients of 238 sediment concentration, e.g. as in the formation of natural lev-239 ees 53 . We find the equilibrium pond size scales with the size of 240 the local basin and increases with the rate R of RSLR ( Fig. 5A 241 left, see Experimental Procedures for pond size calculation). 242 However, isolated ponds tend to connect to the channel net-243 work via the formation of new small channels 41, 42, 49 , thereby 244 increasing channel density and shrinking tidal basins. Based 245 on this, we assume in our model that once ponds are deep 246 enough they connect to channels and become a source of 247 sediment and tidal flow (see Experimental Procedures). Re-248 gardless of the specific conditions for when and how ponds 249 connect, simulations show there is no marsh equilibrium as 250 long as permanent ponds are able to connect to the channel 251 network. Instead, marshes experience a continuous (runaway) 252 fragmentation at a rate controlled by the ratio R/R c (Fig. 5A, 253 right). 254 The runaway fragmentation can be understood as follow: 255 although there are more channels (and connected ponds) 256 to potentially redistribute sediments into the marsh platform, 257 the sediment will be deposited closer to the banks as water 258 flow slows down in the now smaller basins (see Eq. 2). As 259 a result, the drowning threshold R c = A c o + A c i (L) is crossed 260 around the watershed divide of the new system, leading to 261 marsh drowning at ever smaller scales. Therefore, with time, 262 marsh fragmentation propagates from large to small scales 263 following the adjustment of the channel network and tidal flows, 264 until most of the marsh is lost. 265 We can obtain an upper-bound for the threshold rate of 266 RSLR for the onset of runaway marsh fragmentation (R c = 267 A c o + A c i (L), Fig. 6) using a theoretical estimation of the max-268 imum contribution of organic accretion for salt marshes 1 269 (A c o ⇡ 3mm/yr). This value is consistent with accretion rate 270 data of Mid-Atlantic US salt marshes and falls within a broader 271 range of direct and indirect estimations of organic accretion 272 rates of marshes elsewhere (see Fig. S5 and Supplemental 273 Experimental Procedures). Similarly to the trend of inorganic 274 accretion rates with tidal range (Fig. 4), the predicted threshold 275 R c (Fig. 6) shows a fundamental vulnerability for microtidal 276 marshes (d z < 1.5m) and marshes with relatively low sediment 277 supply (average SSC at the channel bank or marsh edge in 278 the range C 0 < 20 g/m 3 ). particularly vulnerable (Fig. 6). 343 We thus provide a mechanistic explanation for the widely 344 observed fragility of microtidal marshes 11 and show this vul-345 nerability is intrinsic and tied to the dominant role of organic 346 accretion. Therefore, factors altering biomass productivity and 347 decomposition, such as eutrophication, elevated CO 2 and cli-348 mate warming 10,11,19,59 , could decide the mid-term response 349 of global microtidal marshes, while measures aimed at increas-350 ing sediment delivery could have limited success. 351 Runaway marsh fragmentation 352 The runaway marsh fragmentation induced by the approximate 353 scale invariance of sediment deposition 44 , constitutes a new 354 form of marsh destabilization that transforms the local cross-355 ing of the marsh drowning threshold into the onset of eventual 356 widespread marsh loss. This mechanism only requires that 357 connected ponds decrease the size of local drainage basins, 358 regardless of whether they deliver sediment to the marsh plat-359 form or not. In the best case scenario depicted in Fig. 5A, 360 connected ponds redistribute inorganic sediment as effective 361 as large channels or mud flats, which is not the case in reality. 362 Any decrease in sediment delivered by connected ponds leads 363 to lower inorganic accretion rates on the surrounding marshes, 364 thereby accelerating marsh drowning. 365 The scale invariance of sediment deposition, where sedi-366 ment is deposited closer to the banks in smaller basins, under-367 pinning the runaway marsh fragmentation is consistent with 368 observations that an increased density of artificial channels 369 does not increase overall sedimentation (e.g. Louisiana 60 ) and 370 in some cases resulted in subsidence (e.g. New England 61 ). 371 Furthermore, the predicted acceleration of marsh fragmen-372 tation with the rate of RSLR (Fig. 5A) is consistent with the 373 rapidly increased rate of historic marsh loss measured in the 374 Mississippi Delta as RSLR accelerated 58 . 375 The marsh fragmentation mechanism explains the formation 376 of a broad range of pond sizes, and predicts that their size 377 distribution should follow a power-law, in agreement with data 378 from Blackwater marshes ( Fig 5B). It also predicts a particular 379 temporal sequence of marsh fragmentation, as large initial 380 ponds eventually lead to smaller ones at a rate increasing with 381 the rate of RSLR relative to the drowning threshold ( Fig. 5A), 382 and suggests the area of the larger ponds depends on the 383 initial distribution of tidal basin areas. This multi-scale mecha-384 nism complements existing models of pond growth driven by 385 lateral expansion instead of RSLR 40, 62 .

387
We derive a simplified model of sediment transport in the 388 absence of erosion that explains patterns of sediment depo-389 sition and marsh vulnerability in a wide variety of conditions. 390 Our model leads to an analytical prediction of inorganic ac-391 cretion that complements direct measurements of accretion, 392 which necessarily reflect historical rather than future environ-393 mental conditions 2 . We predict a new form of marsh destabi-394 lization characterized by a progressive fragmentation of the 395 marsh platform, triggered by the drowning of interior marshes. 396 The threshold for this runaway marsh fragmentation is much 397 lower than existing predictions 13,63 and is largely decoupled 398 from inorganic sediment supply in microtidal environments, 399 which explains the observed fragility of microtidal marshes. 400 Beyond microtidal marshes, the much-lower marsh fragmenta-401 tion thresholds predicted by our model suggest a re-evaluation 402 of the resiliency of global marshes under current and future 403 where the bar denotes an average of the form where t(D) is the rescaled local inundation time and D(x) = 452 d z/2 Z(x) is the local depth. 453 Because the main effect of a non-flat marsh platform is to 454 change the local inundation time t(D), this averaging removes, 455 in a first approximation, the dependence on marsh elevation 456 and thus its solution has the form C ⇡ C(x). Therefore, we can 457 use the numerical solution of Eq. 5 for a flat marsh to obtain a 458 correlation between the average sediment flux per unit width 459 (QC) and the average suspended sediment concentration (C). 460 This correlation is expected when transport is dominated by 461 advection instead of diffusion. 462 Indeed, in the range x/L . 0.6, we find (see Fig. S2) Substituting the advection approximation (Eq. 8) into Eq. 6, 472 we get an equation for the average suspended sediment con-473 centration which has the exponentially decaying solution with decay length From Eq. 6, the scaling of the decay length has the more 478 general form L c µ Q/w f (as can be verified using Q µ d zL/T ), 479 which is equivalent to the scaling of the decay or deposition 480 length in unidirectional turbulent suspensions 50 : L c µ HU/w f µ 481 Q/w f , where H is the flow depth, U is the (constant) flow 482 velocity and Q µ UH is the water discharge per unit width. 483 Finally, the boundary condition C(0) in Eq. 10 is obtained 484 numerically from Eq. 5 by averaging C(0,t) over one tidal cycle, 485 which gives (see Fig. S3) with fitting function This function quantifies the average sediment concentration 488 of the ebb flow leaving the marsh platform. Defining C(0) ⌘ 489 [C flood (0) +C ebb (0)]/2, substituting Eqs. 11 and 12, and using 490 our assumption of a constant concentration at the marsh edge 491 during flood (C flood (0) = C 0 ), we get, For small tidal ranges, the rescaled falling velocity diverges, 493 C ebb (0) ! 0 and most of the sediment is deposited on the 494 marsh. For large tidal ranges, the opposite is true, w + f ! 0 and 495 C ebb (0) ! C 0 , i.e. most of the sediment leaves the march.
In  The inorganic accretion rate A i (x, Z,t) is given by Eq. 14 541 and can be written in terms of the critical accretion rate in the 542 marsh interior, A c i (L) = A i (L, D c ), as: Using Z c /d z = 0.15 as the critical elevation for marshes (corre-553 sponding to D c = 0.35d z, see Fig. 1) we get t(Z c ) = 0.35.

554
The function`(x,t) in Eq. 16 generalizes the concept of the 555 distance x to the marsh edge to account for the formation of 556 new connected ponds. We assume that connected ponds 557 change the geometry of the drainage basin and become a 558 new source of both tidal water and inorganic sediment with 559 concentration C 0 . As ponds get deeper than Z t and connect to 560 the channel network, we update the term`(x,t) to reflect the 561 positions x j of the new marsh edges (defined by the condition 562 Z(x j ) = Z t ), and corresponding watershed divides (defined 563 as the midpoint between neighboring channels or connected 564 ponds.)  575 We use a rescaled critical elevation Z c /d z = 0.15 consistent 576 with field data (Fig. 1B), and assume ponds with a depth 577 around MSL connect to channels, thus Z t /d z = 0. We change 578 the rescaled RSLR rates R/R c in the range 0.8-5. The initial 579 condition is a marsh platform of rescaled elevation Z/d z = 0.4 580 and unit rescaled length, limited by tidal channels at both sides. 581 For the pond size distributions shown in Fig. 5B, we choose a 582 10km domain size. The scale invariance of spatial sediment deposition patterns 585 leads to a similar scale invariance in the size, or diameter L p , 586 of the resulting ponds. Assuming the edge of the pond, a 587 distance x p = L L p /2 from the channel bank, is at equilib-588 rium with RSLR at the critical depth D c , then R = A c o + A c i (x p ) 589 (Eq. 15). Substituting Eq. 16 with Z(x p ) = Z c and rescaled po-590 sition of the pond edge`(x p ) = x p /L = 1 L p /(2L), and using 591 the definition of the drowning threshold R c = A c o + A c i (L), the 592 rescaled equilibrium pond size is where,`c = b /w + f = b d z/(Tw f ) is the rescaled sediment con-594 centration decay length. 595 The rescaled equilibrium pond size (Eq.18) has two limits: 596 no permanent ponds (L p = 0) for R  R c , and no marshes 597 (L p = 2L) above the highest drowning threshold at marsh edge, 598 Fig. 5A). Note that 599

6/15
this pond size is a minimum value as we assume no lateral 600 pond erosion besides marsh drowning.

602
To only test the dependence on the distance to channel, 603 reported accretion rates A i for Phillips Creek ( Table S1 for details). compared to measurements of averaged sediment concentration C (A 34 , B 35 and C 51 ) and inorganic accretion rate A i (D 34 , E 21 and F 52 ) (symbols). A ⇤ i is the depth-corrected accretion rate (see Experimental Procedures for more information). The scaling of the decay length is obtained from the model as L c = 1.5Ld z/(Tw f ) (e.g. Eq. 2), where d z is the tidal range, T is tidal period and w f is the effective sediment falling velocity. In all cases L is taken as the maximum distance to a channel reported in the data, d z (d z ⇤ ) is the reported tidal range (average/typical tidal range during the measurement period), and we use the generic value w f = 10 4 m/s 3, 22, 64 unless stated otherwise. Values of C(0) and A i (0) were fitted to data. Mass accretion rate data was converted to volume accretion rate using an effective density of inorganic sediments deposited in the marsh r i ⇡ 2g/cm 31

13/15
Rate of Relative Sea Level Rise ( )   For R < R c , shallow ponds can recover (bottom center) and marshes reach a non-flat equilibrium state. For R > R c , marsh drowns and form ponds. If those ponds remain isolated, the marsh eventually reaches equilibrium. Otherwise, a self-similar mechanism of pond formation and basin reduction leads to a runaway marsh fragmentation. (B) Exceedance probability distribution of pond areas in Blackwater, MD (representing ponds larger than 50m 2 within the white region in (C), see 55 for details on data acquisition, data available in Table S2)   ) as function of tidal range, for different values of the average suspended sediment concentration at the channel bank C 0 representing typical low, mid and high sediment supply conditions (see Fig. 4B). We use w f = 10 4 m/s and r i = 2 g/cm 3 for the calculation of the critical inorganic accretion A c i (L) (as in Fig. 4), and assume an organic accretion rate A c o = 3mm/yr, consistent with a meta-analysis of field data ( Fig. S5 and Supplemental Experimental Procedures). Symbols represent predictions for specific locations including Blackwater, MD; Plum Island, MA; Phillips Creek, VA and Georgia (we use values shown in Fig. 4B). Current RSLR rates for those locations are in the range 3.5 ± 1.5mm/yr (red line and shaded area). Organic accretion rates in salt marshes are in the range 3.0 ± 2.0mm/yr (green line, shaded area and error-bars, see

k t s G W 4 E d h O F N A o E d o J J Y + 5 3 n l B p H s t H M 0 3 Q j + h I 8 p
A z a q z U a l w 1 B u 6 g X H G r 7 g J k n X g 5 q U C O 5 q D 8 1 R / G L I 1 Q G i a o 1 j 3 P T Y y f U W U 4 E z g r 9 V O N C W U T O s K e p Z J G q P 1 s c e y M X F h l S M J Y 2 Z K G L N T f E x m N t J 5 G g e 2 M q B n r V W 8 u / u f 1 U h P e + R m X S W p Q s u W i M B X E x G T + O R l y h c y I q S W U K W 5 v J W x M F W X G 5 l O y I X i r L 6 + T d q 3 q X V d r D z e V O s n j K M I Z n M M l e H A L d b i H J r S A A Y d n e I U 3 R z o v z r v z s W w t O P n M K f y B 8 / k D q L W N 0 g = = < / l a t e x i t > z = 1m < l a t e x i t s h a 1 _ b a s e 6 4 = " c 2 x N 0 g G i R H s l w C R s H I Z X Q 0 d 9 g l Y = " > A A A B / H i c d V D L S s N A F J 3 4 r P U V 7 d L N Y B F c h a Q W a h d C w Y 3 L C v Y B T S i T y a Q d O p O E m Y k Q Q / 0 V N y 4 U c e u H u P N v n K Y R V P T A h c M 5 9 3 L v P X 7 C q F S 2 / W G s r K 6 t b 2 x W t q r b O 7 t 7 + + b B Y V / G q c C k h 2 M W i 6 G P J G E 0 I j 1 F F S P D R B D E f U Y G / u x y 4 Q 9 u i Z A 0 j m 5 U l h C P o 0 l E Q 4 q R 0 t L Y r L k B Y Q r B O 3 g B n d y V I e T z s V m 3 r X Y B u C S t Z k n a D n Q s u 0 A d l O i O z X c 3 i H H K S a Q w Q 1 K O H D t R X o 6 E o p i R e d V N J U k Q n q E J G W k a I U 6 k l x f H z + G J V g I Y x k J X p G C h f p / I E Z c y 4 7 7 u 5 E h N 5 W 9 v I f 7 l j V I V n n s 5 j Z J U k Q g v F 4 U p g y q G i y R g Q A X B i m W a I C y o v h X i K R I I K 5 1 X V Y f w 9 S n 8 n / Q b l n N m N a 6 b 9 Q 4 s 4 6 i A I 3 A M T o E D W q A D r k A X 9 A A G G X g A T + D Z u D c e j R f j d d m 6 Y p Q z N f A D x t s n / j 6 U R g = = < / l a t e x i t > z = 10m < l a t e x i t s h a 1 _ b a s e 6 4 = " / m + + X Y / L x O J 0 s h t b h E 3 b O c o v 8 J 8 = " > A A A B / X i c d V D L S s N A F J 3 U V 6 2 v + N i 5 G S y C q 5 D U Q u 1 C K L h x W c E + o A l l M p m 0 Q y e T M D M R 2 l D 8 F T c u F H H r f 7 j z b 5 y m E V T 0 w I X D O f d y 7 z 1 + w q h U t v 1 h l F Z W 1 9 Y 3 y p u V r e 2 d 3 T 1 z / 6 A r 4 1 R g 0 s E x i 0 X f R 5 I w y k l H U c V I P x E E R T 4 j P X 9 y t f B 7 d 0 R I G v N b N U 2 I F 6 E R p y H F S G l p a B 6 5 A W E K w R m 8 h I 6 d u T K E 0 X x o V m 2 r m Q M u S a N e k K Y D H c v O U Q U F 2 k P z 3 Q 1 i n E a E K 8 y Q l A P H T p S X I a E o Z m R e c V N J E o Q n a E Q G m n I U E e l l + f V z e K q V A I a x 0 M U V z N X v E x m K p J x G v u 6 M k B r L 3 9 5 C / M s b p C q 8 8 D L K k 1 Q R j p e L w p R B F c N F F D C g g m D F p p o g L K i + F e I x E g g r H V h F h / D 1 K f y f d G u W c 2 7 V b u r V F i z i K I N j c  where both C + (0) and C + (L) are function of tidal range via the rescaled falling velocity w + f (See Fig. S3). Simulation data is shown only for the critical elevation Z c = 0.15d z, but a similar result is obtained for any other elevation.