W HEN DOES FAULTING-INDUCED SUBSIDENCE DRIVE DISTRIBUTARY NETWORK REORGANIZATION ?

Deltas exhibit spatially and temporally variable subsidence, including vertical displacement due to movement along fault planes. Faulting-induced subsidence perturbs delta-surface gradients, potentially causing distributary networks to shift sediment dispersal within the landscape. Sediment dispersal restricted to part of the landscape could hinder billion-dollar investments aiming to restore delta land, making faulting-induced subsidence a significant, yet unconstrained hazard to these projects. In this study, we modeled a range of displacement events in disparate deltaic environments, and observe that a channelized connection with the displaced area determines whether a distributary network reorganizes. When this connection exists, the magnitude of distributary network reorganization is predicted by a ratio relating dimensions of faulting-induced subsidence and channel geometry. We use this ratio to extend results to real-world deltas and assess hazards to deltaic-land building projects.


INTRODUCTION
Deltas provide numerous societal and environmental benefits (Chaplin-Kramer et al., 2019;Hoitink et al., 2020), as well as engender cultural and personal identity (e.g., Zeisler-Vralsted, 2019). Anthropogenic landscape modifications, including levees and channel dredging, limit sediment dispersal necessary to raise deltaic land with respect to water level, thereby promoting widespread land loss by inexorable subsidence (e.g., Paola et al., 2011;Hoitink et al., 2020). Efforts to restore deltaic land (e.g., Peyronnin et al., 2013) rely on a robust understanding of natural sediment dispersal processes (Giosan, Constantinescu, Filip, & Deng, 2013;Temmerman & Kirwan, 2015;Hoitink et al., 2020), as well as how dispersal is affected by disturbances, such as spatially and temporally variable subsidence (Kolb & Van Lopik, 1958). Temporally-discrete faulting-induced subsidence manifests at space and time scales similar to delta restoration efforts (10 3 m and 10 2 yr) Kim, Mohrig, Twilley, Paola, & Parker, 2009), but it is not known whether perturbation by subsidence causes distributary networks to reorganize. Reorganization could limit sediment dispersal extent  or cause sediment to bypass areas targeted for land restoration (Kim et al., 2009). We need to understand the impact of faultinginduced subsidence on distributary network evolution to make robust predictions for deltaic land building.
Deltaic land builds where sediment deposition exceeds accommodation space, and is accordingly dependent on processes that disperse sediment and lower the land surface ; importantly though, dispersal is affected by lowering, and lowering is affected by dispersal. For example, continuous tectonic subsidence "steers" channels toward areas of greater accommodation due to lateral (i.e., cross basin) topographic gradients (Bridge & Leeder, 1979;Heller & Paola, 1996;Kim, Sheets, & Paola, 2010;Reitz et al., 2015). Persistent longitudinal topographic gradients (i.e., down basin) limit channel mobility by "locking" sediment conduits with downstream sinks . In an opposite feedback, tectonic subsidence in deltaic environments commonly includes movement along fault planes, caused by gravitational instability from accumulating sediment (e.g., salt tectonics, slumps, growth faults; . In turn, accumulated growth fault displacement may affect channel orientation and location Gasparini, Fischer, Adams, Dawers, & Janoff, 2016).
These aforementioned patterns constitute time-integrated system behavior and do not predict channel evolution and sediment dispersal at timescales pertinent to coastal restoration (10 2 yr). Indeed, channels ordinarily locate irrespective of time-integrated subsidence patterns (Hickson, Sheets, Paola, & Kelberer, 2005;Kim et al., 2010;Straub & Esposito, 2013;Gasparini et al., 2016;. Significant vertical fault displacement 10 0 -10 1 m over less than one year has been observed , but fault activity varies in space and time (Mouslopoulou, Walsh, & Nicol, 2009;Fossen, 2020). As a result, land building studies have not accounted for faulting-induced subsidence influence on distributary network processes (e.g., Army Corps of Engineers, 2021).
Predicting displacement timing across decades is unlikely (e.g., Geller, 1997), and so site-specific modeling of channel network evolution due to displacement has limited application in coastal restoration. Instead, a framework describing how sediment dispersal and land building are affected if faulting-induced subsidence occurs, would be more helpful. This framework should be interpretable for any delta and across scales of faulting.

SCALES OF FAULTING-INDUCED SUBSIDENCE AND CHANNEL NETWORK REORGANIZATION
Two motivating examples of faulting-induced subsidence on deltas illustrate scales to be considered. At one end, earthquakerelated faulting along ∼10 1 km faults infrequently (∼500 yr period) displaces a moderate portion of a delta area by several meters vertically . At the other end, regional (10 2 km length) and local (10 0 -10 1 km) listric normal faults (growth faults) are ubiquitous on low-lying deltas, and are sporadically but frequently active over centuries, averaging to <1-70 cm/yr . Deltaic faulting-induced subsidence occurs on a spectrum, with bounds on space and time dimensions represented by these scenarios.

Example: large disturbance and wholesale reorganization
The Selenga River delta historical shoreline position documents distributary network reorganization following faultinginduced subsidence. Situated along the margin of rift-basin bound Lake Baikal (Russia), the Selenga River delta is tectonically active (Figure 1a; ). An earthquake in 1862 induced up to 9 m of vertical displacement along the eastern delta edge , lowering >230 km 2 of land below the water surface and creating Proval Bay (Figure 1a; Vologina et al., 2010;. As a result, the shortest path from delta apex to lake level was available via Proval Bay (i.e., eastern deltaic lobe). A bathymetric survey of Proval Bay following the earthquake indicated 1-5 m water depth centrally, and shallower 0.7-1.5 m depth proximal to the delta . In the ∼160 years since, the shoreline along the eastern lobe (i.e., at Proval Bay) has prograded faster than elsewhere along the delta front, and has reached a condition that restores delta planform radial symmetry (Figure 1a; Dong et al., in prep. (Dong et al., in prep.). Arc segments denote the approximate extent of channels within modern deltaic lobes, and labels indicate water flux to each lobe (i.e., a proxy for sediment flux; . b) Mississippi River delta, imaged from Sentinel 2 in March 2021. Faults with down-tosouth motion locate across the delta topset Culpepper et al., 2019), and a subset of these faults reach the surface (inset). Water and sediment exiting main-channel diversions (e.g., at Ironton) may be affected by faulting-induced displacement. Unit-cell water discharge normalized to the inlet flow is overlayed as white-to-red colors (with transparency <5%). Selected simulations from c-f) experiment Set 1 and g-j) Set 2, depicted immediately following displacement event, wherein black contour demarcates fault block extent. Bottom panels (e & f, i & j) demonstrate replicate simulations, where displacement length and block area was unchanged but location varied.
Measurements of Proval Bay water depth and delta subaerial change suggest 6.8 ± 2.1 × 10 7 m 3 of sediment has been delivered to this lobe in ∼160 years (Supporting Information). Scaling deposited volume to the estimated annual bed material input to the delta (2.4 ± 0.2 × 10 5 m 3 ; Dong et al., in prep.) and applying present deltaic lobe partitioning (three lobes, 23% to eastern lobe; Il 'icheva, Pavlov, & Korytny, 2014;Dong et al., 2020, in prep.) or equal lobe partitioning (33% to eastern lobe) accounts for only 10-38% of the observed deposited sediment volume. To match the observed sediment volume, >95% of sediment input to the delta over ∼160 years must have gone to the eastern lobe (Supporting Information), implying wholesale distributary network reorganization following displacement, and subsequently balanced sediment partitioning among deltaic lobes. The 1862 Selenga River delta displacement event exemplifies one end-member scale of faulting, wherein a large disturbance caused a dramatic response in the delta network.

SUBSIDENCE ACROSS A RANGE OF SCALES
Inspired by the compelling example of faulting-induced distributary reorganization on the Selenga River delta and a need to evaluate hazards to Mississippi River delta sediment diversions, we ask: what scale of displacement (vertical slip and areal extent) drives distributary network reorganization, and are equilibrium sediment distribution processes eventually restored? We implemented numerical experiments as a sensitivity analysis to address these questions; we apply a perturbation with varied vertical displacement length (H σ ) and fault-block area (A σ ) to a consistent initial state, and observe subsequent deltaic system evolution. Simulations used the DeltaRCM numerical delta model , implemented in Python as pyDeltaRCM (Moodie, Hariharan, Barefoot, & Passalacqua, in review). DeltaRCM has been robustly validated (Liang, Ge-leynse, Edmonds, & Passalacqua, 2015b;) and used to examine delta morphology and evolution under various external forcings, including sea-level rise , and presence of vegetation (Lauzon & Murray, 2018), and ice and permafrost (Lauzon, Piliouras, & Rowland, 2019;Piliouras, Lauzon, & Rowland, 2021). The perturbation was modeled as instantaneous vertical land movement (i.e., displacement without hanging-wall rotation or translation), with channel evolution and geometry before and following displacement emerging based on model boundary conditions (Supporting Information). We included replicate simulations for each displacement size (after  and simulations without displacement, to assess autogenic variability and quantify uncertainty in distributary response. Experiment Set 1 begins with a deltaic system qualitatively similar to the Selenga River delta (Figure 1a and 2a; parameters listed in Supporting Information). Namely, the model delta maintains several active distributary channels (6-8;Il'icheva et al., 2014; and an approximately radially symmetric subaerial planform (Reitz & Jerolmack, 2012;Dong et al., in prep.). Set 1 varied H σ = {0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5} meters while holding fault block width (B σ ) and length (L σ ) fixed at 6.5 and 12 km, respectively, and locating the fault (i.e., block edge nearest the delta inlet) randomly along an arc 6 km from the delta inlet (Figure 2cf; seven replicates, total 70 simulations). Set 1 displacement lengths and block width range Selenga River delta observations Dong et al., in prep.).
The deltaic system beginning Experiment Set 2 is qualitatively similar to the Mississippi River delta (Figure 1b and 2b; parameters listed in Supporting Information). The Set 2 delta system maintains 1-2 active distributary channels, and an asymmetric planform with high topography isolated along relict alluvial ridges and surrounded by shallow embayments. Set 2 varied displacement length H σ = {0.1, 1, 5} meters and fault block area, by modulating block width B σ = {3.6, 12, 24} km and holding block length L σ = 9.6 km. The fault location was constrained to 30-56 km from the delta inlet and placed randomly ( Figure 2c-i). Subsequent planform change indirectly records distributary network response, whereas network response is directly recorded in spatiotemporal sediment flux partitioning across the landscape. Flux partitioning measurements are rare from deltas , but accessible in simulations. Therefore, we analyzed model simulations with indirect and direct metrics (Supporting Information), ensuring our analyses are comparable to historical records and can be generalized to predict reorganization in real-world deltas. ). Values less than unity describe asymmetry due to faulting-induced subsidence in simulations. c) Sediment flux asymmetry (R Q s ; Equation 4). Time-averaged sediment partitioning is balanced ≈ 0.33, and shift towards unity indicates sediment partitioning due to network reorganization. In a-c), solid lines and error bars are mean and standard deviation for replicate simulations from each displacement length, gray shaded region is autogenic variability without faulting-induced subsidence, and black vertical dashed line marks displacement timing.

Planform morphology
Keeping in mind the objective to build new land by sediment diversion, we quantified subaerial delta land area (A L ) as: where i and j are model cell indices, L and W are the number of cells in each dimension of the model domain, and where dx is the model grid cell size, η is the bed elevation, and H W L is the basin water level. Subaerial delta area abruptly decreases proportionate to displacement length in Experiment Set 1 (Figure 3a), but area does not change predictably for Set 2 (i.e., within autogenic variability; Figure 4a). Deltaic land area unsteadily increases in Set 1 and Set 2 simulations, though Set 2 increase is non-monotonic and more variable (Figures 3a and  4a).
Approximately one-third of the Selenga River delta width was lowered by the 1862 earthquake, creating a radially asymmetric delta planform Dong et al., in prep.). We quantified this delta asymmetry (R r ) in simulations as a ratio: where r σ and r are the mean distance from delta apex to points along the shoreline within the subsided block and across the entire delta, respectively (shoreline identified by the Opening Angle method; Shaw, Wolinsky, . Set 1 deltas become asymmetric immediately following displacement (R r < 1 and outside autogenic variability) for displacement greater than ∼0.2 m ( Figure 3b). Then, R r increases, at a different rate for each displacement length, until the shoreline position reaches a dynamic equilibrium (i.e., delta symmetry, R r ≈ 1; Figure 3b). Fault block area and shoreline do not necessarily overlap in Experiment Set 2, rendering delta asymmetry uninformative for this set.

Distributary network organization
Simulations yield sediment transport data at spatiotemporal resolution unparalleled in real-world systems (Hoitink et al., 2020) and invaluable for understanding delta evolution following faulting-induced subsidence. We quantified shifting distributary network organization via sediment flux asymmetry (R Q s , i.e., sediment partitioning): where Q s,σ is sediment flux into the subsided block, and Q s is sediment flux across the delta, calculated through a circular cross-section positioned at the fault (Figure 2c-i). In Experiment Set 1, sediment flux partitioning is approximately symmetric for displacement lengths ≤ 0.1 m (R Q s ≈ 0.33, as the subsided block measures one-third of the delta width). For larger displacement lengths, sediment is consistently partitioned into the subsided block (R Q s → 1; Figure 3c), and duration of asymmetric flux scales with displacement length (Figure 3c) and is coincident with delta asymmetry (Figure 3b). For Set 2, flux asymmetry following the displacement event depends weakly on displacement length, but decorrelates a short time later (Figure 4b). Shifted flux in Set 1 and Set 2 demonstrates that disequilibrium in deltasurface gradients created by faulting-induced subsidence drives distributary network reorganization to restore equilibrium.
In many Set 2 simulations, fault location is disconnected from the channel network (Figure 2g-j), which categorizes simulations as either "connected" or "not connected", based on whether the subsided block had > 5% of the inlet unit-cell discharge at the time of displacement (Figure 2). Connectedness dictates when sediment first enters the subsided block (T 0 , Figure 4c), with blocks receiving sediment a median time of 9 and 381 years later for connected and not-connected locations, respectively (p = 7.0 × 10 −13 for an unequal variances t-test). Volumetric sediment flux into the subsided block at time T 0 (Q s,0 , Figure 4d) is indistinguishable between connected and not-connected Set 2 simulations (p = 2.3 × 10 −1 ). All Set 1 simulations are connected (> 5% inlet unit-cell discharge), causing networks to swiftly reorganize after displacement (i.e., T 0 ≈ 0), and correlating volumetric sediment flux to the subsided block (i.e., Q s,0 ) with displacement length (Figure 3c).

Style and timing of distributary network reorganization
Differences in reorganization predictability between Set 1 and 2 are due to distributary network configurations unique to each set, that arise from autogenic sediment dispersal characteristics. Set 1 simulations shift sediment dispersal over time via soft avulsion (Edmonds, Paola, Hoyal, & Sheets, 2011), whereas sediment dispersal shifts in Set 2 simulations by deltaic lobe-switching avulsion (Supporting Information;Slingerland & Smith, 2004). Soft avulsion leads to multiple active distributary outlets in a branching channel network , which ensures connectedness (and thus reorganization). In contrast, lobe-switching avulsion abandons channel courses and limits active sediment transport extent, such that connectedness (and thus reorganization) depends on fault location.
Avulsion behavior continues to control spatiotemporal sediment distribution after initial reorganization (or stasis). For example, sediment flux timing of not-connected simulations (T 0 ) is distributed approximately uniformly in time (Figure 4c), which is consistent with periodic avulsions (e.g., Mohrig, Heller, Paola, & Lyons, 2000), and channel location is determined by pre-existing topography (e.g., Hajek & Edmonds, 2014) rather than being steered towards the subsided region. Moreover, indistinguishable sediment flux Q s,0 (Figure 4d) indicates avulsed channels are self-formed (Slingerland & Smith, 2004;Hajek & Edmonds, 2014), because channels convey consistent sediment flux, regardless of displacement scale or timing. Gradual recovery of balanced sediment partitioning following abrupt reorganization in Experiment Set 1 (Figure 3c distributary channels and soft avulsion are observed on the Selenga River delta (similar to Set 1; Dong et al., 2016, in prep.). Additionally, lobe-switching avulsion commonly controls sediment distribution in single-threaded deltaic systems like the Mississippi River delta (similar to Set 2; Roberts, 1997;Jerolmack & Mohrig, 2007). Ephemeral asymmetric shoreline progradation focused along the subsided block in Set 1 (Figure 3c) is consistent with the Selenga River delta historical record (Figure 1a) and modern delta partitioning . Behavioral similarity between model and real-world across disparate deltaic systems is encouraging for efforts to extrapolate model results to inform land building and distributary networks in real-world deltaic environments.

Impacts on land building
Land building at the delta scale is not significantly impacted by faulting-induced subsidence (Figures 3a and 4a). Bolstered by Mississippi River delta sediment budgets indicating sufficient resources for land building (Kim et al., 2009;Nittrouer & Viparelli, 2014;Sanks et al., 2020), we anticipate that faultinginduced subsidence will not stymie delta-wide land building projects there, though land emergence may be delayed (e.g., Figure 3a). Additionally, sediment builds land while focused by reorganization and channels are mobile after filling accommodation space (Figure 3a-b), ensuring that sediment does not bypass the delta topset indefinitely (Figures 3 and 4; Supporting Information). Over time, these sediment dispersal processes (i.e., avulsions) leads to delta-scale landscape sustainability (Passalacqua, Giosan, Goodbred, & Overeem, 2021).
Faulting-induced subsidence and restricted sediment dispersal potentially stymie land building downstream of diversions (i.e., at a local scale; Figures 3c and 4b; Supporting Information). Land-building models not accounting for reorganization (e.g., Army Corps of Engineers, 2021) potentially overestimate extent and location for land building. Though, reorganization is ephemeral, which is favorable for diversion planning; distributive sediment dispersal is restored without intervention ( Figures  3c and 4c-d). Faulting-induced subsidence is a hazard to land building at scales relevant to coastal restoration, but more research is needed to resolve impacts in detail or at specific sites.

Impacts on distributary network organization
Association between paleo-channels and faults mapped in the subsurface led  to hypothesize that channel reorganization depends on the relative scale of fault displacement and channel size. Building on this idea, we determined a dimensionless number relating displacement event and distributary network characteristics that we term the displacement magnitude ( σ ): where F w is fraction of delta width channelized (akin to nearest edge distance; Edmonds et al., 2011), B c is channel width, and H c is channel depth, all measured at the fault distance from delta apex (Supporting Information). Figure 5 shows the relationships between displacement magnitude ( σ ) and distributary network reorganization (∆R Q s ), which we quantified as the difference between sediment flux asymmetry before and after displacement ( Figures 3c and 4b); not-connected simulations from Set 2 were omitted.
Expected reorganization increases as faulting scales larger relative to the distributary network ( Figure 5). Equation 5 codifies variables needed to characterize the distributary system (network density, channel width, channel depth), and so enables generalizing results, with displacement magnitude estimates for real-world distributary systems ( Figure 5, Supporting Information). (Figure 5). For example, the Selenga River delta has a σ range indicative of substantial distributary network reorganization following large and periodic earthquake-related displacements, and the same should be expected for similarly scaled systems. At these scales, reorganization may be recorded in the chemical and sedimentary stratigraphy of the delta (Dong et al., 2016) as unconformities and hierarchical packages that appear similar to autogenic avulsions (Ganti, Lamb, & Chadwick, 2019;Ganti, Hajek, Leary, Straub, & Paola, 2020). At the opposite spectrum end, the Mississippi River delta main-channel σ range indicates minor reorganization ( Figure 5). Intuitively, reorganizing Mississippi-sized channels would require larger displacement than has been documented on the delta .

End-member deltaic configuration examples motivating this study bound scales of distributary network reorganization
At the scale of planned sediment diversions (i.e., scaled to Cubit's Gap), faulting-induced subsidence would drive distributary reorganization and restrict sediment dispersal ( Figure 5).
Note that sediment flux reorganization documented herein is different from channel "locking" , because displacement and channel mobility are not coeval, and channels are mobile after accommodation has been filled. Cubit's Gap scaling informs land building efforts on the Mississippi delta, and shows how the relationship in Figure 5 can be used as a tool to understand potential diversion hazards on deltas globally. At these intermediate σ scales, the distributary network response to faulting-induced subsidence is uncertain ( Figure 5) and will depend on channel location with respect to faulting (Figure 4c-d). As a result, an ensemble modeling approach is necessary to understand the probable range of land building scenarios.
Validating models will require field-based flux partitioning observations from deltas . These data will support an understanding of delta response to perturbations generally, and promote sustainable landscape management (Passalacqua et al., 2021). In this context, we suggest targeting mapped fault locations with temporally-resolved morphodynamic and hydrodynamic measurements, thus providing a baseline for future faulting-induced subsidence scenarios. These data will also inform how diffusive processes such as waves and tides work to fill subaqueous topographic lows without channelized flow (e.g., Obelcz, Xu, Bentley, O'Connor, & Miner, 2018).

DATA AVAILABILITY
This study utilized open-source projects pyDeltaRCM (v1.2) and DeltaMetrics (v0.1). The project workflow is available on GitHub at https://github.com/amoodie/paper _resources under folder "Moodie_faultingsubsidence". The GitHub repository includes scripts to configure and execute model runs, process and analyze data, and generate figures in this article. Additional data will be made available in a public repository if paper is accepted. 1. Additional methods description 1.1. Mass balance into Selenga River delta eastern lobe following 1862 earthquake A compelling example of distributary network reorganization due to faulting-induced subsidence comes from the Selenga River delta (Russia). The Selenga River delta is situated along the margin of Lake Baikal, a lake bounded by tectonically active intracontinental rift-basin faults (Figure 1a; . The modern delta channel network divides into three distinct deltaic lobes (Figure 1a; Ilicheva et al., 2014), that each receive a moderate proportion of the water and sediment entering the delta (25-50% each; Figure  1a; Dong, 2020).
An earthquake in 1862 induced up to 9 m of vertical land subsidence along the eastern delta edge , lowering more than 230 km 2 of land below the water surface and creating Proval Bay (Figure 1a; Vologina et al., 2010;; lowered land comprised unconsolidated Quaternary and older Tsagan Steppe deposits . As a result, the shortest path from delta apex to lake level was available via Proval Bay (i.e., the eastern delta lobe). A bathymetric survey of Proval Bay following the earthquake indicated 1-5 m water depth centrally, and shallower 0.7-1.5 m depth proximal to the delta front . In the ∼160 years since, the shoreline along the eastern lobe (i.e., at Proval Bay) has prograded faster than elsewhere along the delta front (Dong, 2020), and has reached a condition that restores delta planform radial symmetry.
We are not aware of historical maps providing adequate temporal resolution to quantify distributary network reorganization following the earthquake. Instead, we apply a mass balance framework, including uncertainty via Monte Carlo sampling (number of samples n = 1000), to indirectly discern whether reorganization was likely following the earthquake. We use published data and measurements from modern satellite imagery to determine that the Selenga River delta distributary network very likely reorganized, with water and sediment flux focused in the eastern delta lobe, following the 1862 earthquake. The following paragraphs elaborate on the approach, and mass-balance input parameter distributions.
The volume of sediment deposited in the ∼160 years since the earthquake was determined by multiplying an estimated deposit area (D A ) by an estimated deposit thickness (D Z ). For D A , we used the areal extent between the modern shoreline and the shoreline as mapped in ∼1862 (Dong, 2020), and determined D A = 67.1 ± 12.0 km 2 (Figure S1a); we assessed areal uncertainty by mapping the polygon at varied resolutions. We then used this polygon to constrain the area over which to estimate deposit thickness. Deposit thickness was determined as the water depth immediately following the earthquake (i.e., ignoring effects of delta or basin slope); measurements of Proval Bay water depth are determined from data published in Vologina et al. (2010) and originally derived from Orlov (1872). These measurements describe depth range in Proval Bay proximal to the delta front 0.7-1.5 m, and the mean and standard deviation used for the mass balance of D Z = 1.1 ± 0.3 m ( Figure S1b). Together, estimated deposit areal extent and thickness yield a volumetric sediment deposit D V = 6.8 ± 2.1 × 10 7 m 3 to the eastern delta lobe in a duration D T =150-160 years; we used D T ∼ U (150, 160) for Monte Carlo samples.
Sediment flux to the delta has likely varied over the last century, and is poorly constrained (e.g., Chalov et al., 2016). We used an estimate of bed-material load from Dong (2020) of Q bm = 2.35 ± 0.22 × 10 5 m 3 /yr for annual sediment input to the delta. For simplicity, we used a fixed depositional porosity φ = 0.25 (Leopold et al., 1964).
We define a governing mass-balance equation for a single deltaic lobe as: where V is deposit volume, T is a time duration, and F is the fraction of sediment input to the delta that is partitioned to the deltaic lobe of interest. Equation S1 can be rearranged and solved with different constraints. In the second row of Figure S1 (d-f), we set V equal to D V and solved for T , varying the sediment partitioning to assess how much time would be required to reproduce the observed sediment volume. We varied the present delta partitioning according to the modern partitioning F = 0.23 (23% to eastern lobe), equal delta partitioning F = 0.33 (33% to eastern lobe), or focused delta partitioning F = 1.0 (100% to eastern lobe). The third row of Figure S1 (g-i) depicts the results discussed in the introduction of the main text, where we solved Equation S1 for V and taking a ratio V /D V , while setting T ∼ D T and varying delta partitioning as F =0.23, 0.33, and 1.0.
Finally, we determined the flux partitioning needed to produce the observed deposit volume by setting T ∼ D T and V = D V and solving for F. From this calculation, we determine that 95-189% of the sediment input to the delta is needed to reproduce the observed record. Therefore, and in conjunction with highly asymmetric shoreline progradation (Dong, 2020), we infer that it is very likely that substantial distributary network reorganization occurred following the 1862 earthquake on the Selenga River delta.

Simulation pyDeltaRCM parameters
Set 1 and Set 2 parameters are highlighted in Table S1. Model parameter sets identified by previous studies informed boundary conditions chosen for Set 1 and Set 2 simulations . To qualitatively match the Selenga River delta with Set 1 simulations, we imposed a relatively high proportion of sand input to the delta (55%)  and maintained a moderate basin depth at 5 m. High sand proportion coupled with fixed sea level produced an approximately radially symmetric delta planform with local embayments around former distributary channels (Figure 2a). In Set 2 simulations configured to qualitatively match the Mississippi River delta, we imposed a low proportion of sand input (35%), slow sea level rise at 2 mm/yr, and a deep basin at 10 m ( Figure 2b). Together, these parameters led to a condition that drowned inter-distributary bays between avulsions ( Figure  2b). Additional parameters (Table S1) scaled the Set 1 and Set 2 deltas and influenced the size of channels.

Graphical depiction of metrics
This Section includes graphical depictions of the metrics (equations) described in described in the main text Section 4 ( Figure S2). Figure S2a shows a simulation from Set 1 (H σ = 2) immediately following the displacement event.
We define delta length asymmetry (R r ) as: (i.e., Equation 3) where r σ and r are the mean Euclidean distance from delta apex to points along the shoreline within the subsided block and across the entire delta, respectively. Figure  S2b depicts the shoreline as determined by the Opening Angle method , setting the shoreline contour threshold to 75 • . Given the set of cells along the shoreline I and a set of cells coincident with the subsided block S, we find the shoreline cells inside the subsided block I S = I ∩ S (darker shoreline segments in Figure S2b) and other shoreline points I X (lighter shoreline segments in Figure S2b). Then, r σ and r are calculated as: where |I S | and |I X | are the length of sets I S and I X , x i and y i are the xand y-coordinate of the i th index in the set, and x 0 and y 0 are the xand y-coordinate of the delta inlet. Figure S2c shows the circular section used for computing the flux asymmetry metric, which we define in the main text Section 4 as: (i.e., Equation 4) where Q s,σ and Q s are the sediment flux into the subsided block and across the delta, respectively. We computed Q s,σ and Q s as: where q s is the unit-cell directed sediment discharge (m 2 /s), and {A, B,C, D} are ordered points defined on the along-section coordinate s of the circular section ( Figure S2c). The definite integral was evaluated numerically via the trapezoidal rule. Finally, the subaerial delta land area (A L ) was quantified as: (i.e., Equation 2) where dx is the model grid cell size, η is the bed elevation, and H W L is the basin water level elevation. Figure S2d shows the area of subaerial land by covering the excluded area with a black shading. Note that channels are usually excluded from the calculated subaerial land, because the bed elevation in these locations is below water level (Equation S9).

Estimating displacement magnitude for real-world deltas
We estimated the displacement magnitude ( σ ) for real-world deltaic systems along a cross section placed approximately halfway between the delta apex and coast ( Figure S3). We acknowledge that this placement is somewhat arbitrary, given that our study identifies the importance of fault location with respect to the channel network. However, we expect that placement in the central delta will yield representative displacement magnitude values.
We relied on published measurements of fault displacement and channel width and depth to constrain displacement magnitude ( σ ) for real-world deltaic systems (Esposito et al., 2013;Nittrouer et al., 2008Nittrouer et al., , 2012Dong et al., 2019;. We selected subsets of data, as possible, to best characterize the channel geometries within the central delta region. We measured the fraction of delta width channelized from Sentinel 2 imagery collected in 2020 and 2021 for the Selenga and Mississippi River deltas, respectively ( Figure S3). Using this imagery, we identified channelized lengths along the cross section, and then simply computed the fractional length channelized.
Channel geometry and network data are compiled in a spreadsheet available within the Github repository https://github.com/amoodie/paper resources under folder "Moodie faultingsubsidence". Shapefiles for measurements of fraction channelized delta width are also included there.

Selected simulation results
In this section we include planform view of the bed morphology and water discharge field following displacement, from selected representative model simulations. Selected simulations from Experiment Set 1 depict the range of behavior in reorganization following displacement, as well as how soft avulsion maintains multiple active distributary channels across the delta topset. Selected simulations from Experiment Set 2 show how reorganization depends on "connectedness" (see main text), and the nature of lobe-switching avulsions in this experimental set. See Figure S4 and S5 captions for detailed description. 0.80 ± 0.35 i) Figure S1. Mass balance framework inputs and results; distributions represent samples from Monte Carlo assessment with n = 1000 samples. Inputs to mass balance framework a) deposit area and b) deposit thickness. d-f) Results from solving Equation S1 for time needed to reproduce the deposit for different sediment partitioning fractions. g-i) Results from solving Equation S1 for fraction of observed deposited volume that would be deposited over historical record duration for different sediment partitioning fractions.

a) b)
A D B C c) d) Figure S2. Graphical depiction of a) planform delta morphology following displacement, and metrics computed in main text (b-d). b) Shoreline identified by the Opening Angle Method , split into segments within the subsided block (dark pink) and elsewhere along the delta (light pink); i.e., depiction of Equation 3. c) Circular cross section along fault for sediment flux asymmetry, where flux to subsided block is computed along section between points B and C (dark pink), and delta-wide flux is computed along section between points A and D (light pink); i.e., depiction of Equation 4. d) Subaerial delta area depicted with areas where the bed elevation is below the water level; i.e., depiction of Equation 2.  Figure S4. Planform evolution of the bed morphology and water discharge field following displacement for selected simulations from Experiment Set 1. In simulation depicted in first column (a-e), displacement length is relatively small and does not induce reorganization. Several soft avulsions occur throughout the 88 yr period; note that multiple distributary outlets remain active at all times. Simulation depicted in second column (f-j) exemplifies the subtle reorganization induced by a moderate displacement length. For ∼44 yr, flow is towards the subsided region, before avulsion relocates the bulk of the water discharge. Simulation depicted in third column (k-o) exemplifies the most extreme case, where displacement length is large, and water discharge is focused in the subsided region for the duration of the simulation. Sediment discharge follows water discharge, and a sub-delta builds within the bay created by displacement; at the end of the simulation, the delta is symmetric.  Figure S5. Planform evolution of the bed morphology and water discharge field following displacement for selected simulations from Experiment Set 2. In simulation depicted in first column (a-e), the fault area is relatively small and "not connected" to the inlet. The distributary network does not reorganize, and when a lobe-switching avulsion occurs later at 45-60 yr, the new channel pathway is not towards the subsided block. Simulation depicted in second column (f-j) exemplifies reorganization due to faulting-induced subsidence. Following displacement, the main distributary channel swiftly redirects to the subsided area, where it remains until the area has been filled with sediment and a lobe-switching avulsion relocates flow (∼45 yr). Simulation depicted in third column (k-o) shows reorganization, but not as clearly or substantially as second column. Difference between third and second column shows how reorganization depends on fault area (these simulations had the same displacement length and were "connected"). At the end of depicted time in third and first columns (60 yr), channel path selection during avulsion is towards the left of the frame, indicating that topography controls the pathway (Reitz et al., 2010).