Land loss in the Mississippi River Delta: Role of subsidence, global sea-level rise, and coupled atmospheric and oceanographic processes

The Mississippi River Delta in coastal Louisiana has suffered large-scale land loss during the historic period and is representative of a global phenomenon where low-elevation deltaic coasts are increasingly at risk because of disrupted sediment supply and accelerated global sea-level rise. Land loss is a natural part of deltaic evolution over time, and most of the land loss in the Mississippi River Delta occurred after individual delta-plain headlands were abandoned as active constructional landscapes, but before 1932 when collection of air photos would make repeat land loss measurements possible. A coastwide land loss of ~5000 km 2 is now well documented for the period 1932 to 2016, which corresponds to a mean rate of ~57 km 2 yr (cid:0) 1 . We use a LiDAR digital topobathymetric model to hindcast land-area changes through time for 1950 – 2010 by incrementally restoring elevation lost due to subsidence, global sea-level rise, and annual anomalies in mean sea level. Our results support the view that the magnitude and spatial distribution of 20th century land loss can be explained by an unfortunate convergence of ongoing subsidence, greatly reduced sediment dispersal due to levee construction, and acceleration of global sea-level rise. Other factors have contributed to land loss on local scales, but the magnitude of land loss that has occurred would have occurred anyway due to subsidence, lack of sediment input, and accelerated sea-level rise. Multidecadal accelerations and decelerations in land loss from 1950 to 2010 have been observed, and attributed to accelerations and decelerations in subsidence due to subsurface fluid withdrawals. However, non-linear land loss represents measurements that were made when water levels varied due to annual to multidecadal anomalies in mean sea level. Anomalies in mean sea level are driven by flux of water into and out of the Gulf of Mexico from the Atlantic, as well as the Atlantic Multidecadal Oscillation (AMO), which produces precipitation anomalies in the Gulf of Mexico drainage area, anomalies in Mississippi River discharge to the Gulf of Mexico, and changes in wind directions that serve to trap water along the coast and elevate coastal sea level, or advect water away from the coast to lower coastal sea level. Sea-level anomalies of the scale described here amplify or suppress the secular trend of global sea-level rise and its impacts on low-elevation delta plains as they respond to ongoing subsidence and anthropogenic disruption of sediment dispersal.


Introduction
The Mississippi River Delta (MRD) in coastal Louisiana (Fig. 1) has suffered large-scale conversion of land to water (hereafter land loss) during the historic period, and is representative of a global-scale phenomenon where low-elevation deltaic coasts are increasingly at risk because of disrupted sediment supply and accelerated global sea-level rise (Syvitski et al., 2005;Syvitski, 2008;Syvitski et al., 2009;Nienhuis et al., 2020).Land-loss in coastal Louisiana is well documented (e. g., Gagliano et al., 1981;Penland et al., 1990;Britsch and Dunbar, 1993;Boesch et al., 1994;Twilley et al., 2016;Couvillion et al., 2017), and the economic and social implications are now widely recognized (e.g., Day Jr et al., 2007;LCPRA (Louisiana Coastal Protection and Restoration Authority), 2017).The most comprehensive and detailed estimates of total land loss are derived from a time series of air photos and satellite imagery that span the period 1932-2016, and show a total land loss 4833 km 2 (Fig. 2a), which corresponds to a century-scale mean rate of ~57 km 2 yr − 1 for the MRD as a whole (Couvillion et al., 2017).
In this paper we further interrogate MRD land loss, including the issue of land loss accelerations and decelerations.We first provide an overview of delta formation and the variables that drive land loss, and place 20th-century land loss within a late Holocene context.We then use regional-scale LiDAR digital topobathymetric data to (a) estimate sensitivity of the MRD to changes in water level, (b) hindcast land loss that can be attributed to sea-level rise and subsidence, and (c) estimate the amount of sediment that would have been required to sustain deltaplain surface area.Last, we propose a chain of causality for multidecadal accelerations and decelerations in land-loss rates that features anomalies in mean sea level driven by annual to multidecadal changes in the coupled ocean-atmosphere system.Although land loss measurements in Couvillion et al. (2017) extend back to 1932, we restrict our analyses to 1950-2010 because of the availability of atmospheric, hydrologic, and oceanographic data needed for our analyses.

The late Holocene to 20th century Mississippi River Delta Region
Large rivers like the Mississippi construct deltaic headlands that extend into the coastal oceans.Active headlands have primary distributary channels that bifurcate downstream to deliver sediment that aggrades the delta plain and progrades the delta front when sediment dispersal to the delta exceeds the accommodation produced by relative sea-level rise (the product of subsidence and sea-level change), then degrade when the reverse is true (Muto and Steel, 1997;Paola et al., 2011).Over millennia, deltaic landscape evolution follows what has been referred to as the delta cycle (Fig. 3) where an active deltaic headland builds seaward, then avulsion relocates the primary distributary channel and river mouth to construct a new headland laterally along the coast.The abandoned delta front is then eroded, and sediment is redistributed by waves and currents, and the abandoned delta plain degrades, subsides, and is ultimately submerged (Gagliano and van Beek, 1975;Coleman, 1988;Penland et al., 1988;Roberts, 1997).
The Holocene MRD extends ~300 km along the northern Gulf of Mexico (GoM) coast and spans an area of >35,000 km 2 .Like most of Earth's major river deltas, the MRD formed over a late Pleistocene glacial-period incised valley that filled with sediment during post-glacial sea-level rise, then developed expansive aggradational and progradational deltaic landscapes as post-glacial global sea-level rise decelerated (Stanley and Warne, 1994).The Holocene MRD consists of 6 episodes of deltaic headland construction that formed successively over the last ~7000 yrs., but the extant Holocene landscape is comprised of the late Holocene St. Bernard (ca.4000-1700 yrs.BP), Lafourche (ca.1700-500 yrs.BP), and Plaquemines-Balize (ca.1400 yrs.BP to present) deltaic headlands, and the embryonic 20th-century to present Atchafalaya-Wax Lake deltas (see Fig. 1) (Frazier, 1967;Törnqvist et al., 1996;Roberts and Coleman, 1996;Roberts, 1997;Chamberlain et al., 2018).The MRD is linked to the Louisiana Chenier Plain, which extends west from the Holocene delta to the Louisiana-Texas border and consists of sandy and shelly ridges of late Holocene age that are separated from each other by mud flats, swamps, and lakes (Gould and McFarlan, 1959;Penland and Suter, 1989;McBride et al., 2007;Hijma et al., 2017).
At a more detailed level, the MRD delta plain includes expansive wetlands that are especially sensitive to the balance between deposition of mineral and organic sediment and the loss of elevation due to subsidence and sea-level rise, a balance formulated in Paola et al. (2011): where: A w = land-surface area of the delta region that can be sustained, Q s = the volumetric rate of sediment supply, fr = the fraction of sediment trapped in the delta topset, r o = the ratio of organic matter to mineral sediment in the sediment column, C 0 = the fraction of solids in the sediment column (1-porosity), σ + H = relative sea-level rise (subsidence rate plus regional sealevel rise) Within this context, wetland sustainability can be framed in terms of elevation capital, defined as marsh elevation relative to the lower threshold for vegetation growth and the intertidal zone (e.g., Reed, 2002;Cahoon et al., 2019).Elevation capital increases over time if deposition rates exceed relative sea-level rise, and decrease over time if the reverse is true: Reed et al. (2020) argue that MRD marshes are resilient to flooding from daily tidal cycles, but flooding of brackish and saline marshes for a period of only 2 years due to sea-level rise is sufficient to result in inundation collapse and conversion of the marsh to water.Moreover, while our research focuses on the vertical component of land loss due to inundation, marshes with low elevation capital can also exhibit significant edge erosion (Fagherazzi et al., 2013;Mariotti, 2020).For example, sea-level rise is thought to trigger tidal channel widening via edge erosion, and may be especially important in settings with high rates of relative sea-level rise, low tidal ranges, and limited sediment supply (Mariotti, 2020).

Late Holocene to 21st century sea-level change
During the Last Glacial Maximum (ca.26-19,000 yrs.BP), global Fig. 3.The delta cycle concept (redrawn and modified from Penland et al., 1988, Blum andRoberts, 2012, andearlier work), illustrating the autogenic evolution of deltaic headlands following avulsion and relocation of feeder distributary channels.(a) Map views of changes through time in morphology and environments as a deltaic headland degrades.(b) Trends in deltaic headland area, shoreline length and other characteristics as deltaic headlands form then decay over time, with the state of individual MRD Holocene deltaic headlands as shown.mean sea level (GMSL) was ~130 m below present, then deglaciation resulted in GMSL rise, with ocean water volumes approaching present levels by ca.7-6000 yrs.BP (e.g., Lambeck et al., 2002;Lambeck et al., 2014;Peltier et al., 2015).Kopp et al. (2016) and Walker et al. (2021) summarize GMSL change from the last ~2000 yrs.(since 0 CE) based on empirical data from >30 localities around the world and conclude the amplitude of GMSL variability was ±70 to ±110 mm prior to the era of instrumentation, with annual rates of change estimated to have been less than ±0.3 mm yr − 1 until the late 19th century.The 20th to 21st century GMSL record has been reconstructed from a global network of tide gauges, and from satellite altimetry since 1993 (see Church and White, 2011;Church et al., 2013;Hay et al., 2015;Cazenave et al., 2018;Nerem et al., 2018;Dangendorf et al., 2017Dangendorf et al., , 2019;;Frederikse et al., 2020).As summarized in Oppenheimer et al. (2019), GMSL rise was ~1.5 ± 0.4 mm yr − 1 from 1900 to 2010, but accelerated by − 0.002-0.019mm yr − 1 over this period to ~3.22 ± 0.37 mm yr − 1 during 1993-2015.Oppenheimer et al. (2019) conclude that 20th and early 21st century rates of GMSL rise are several times higher than during the late Holocene, which is the period during which the MRD St. Bernard and Lafourche deltaic headlands were constructed.including melting of land ice, thermal expansion of the oceans, ocean circulation, and wind and pressure changes, but also the presence or absence of vertical land motion (VLM), especially subsidence.Reconstructions of late Holocene RSL change for the MRD constrain rates of rise to <0.5 mm yr − 1 over the last ~2000 yrs., decelerating until the 1800s (Törnqvist et al., 2004;González and Törnqvist, 2009;Törnqvist et al., 2020).These reconstructions excluded, by design, contributions to RSL rise from subsidence due to compaction of Holocene deltaic sediment (see below), and are interpreted to be dominated by ongoing glacial-isostatic adjustment.20th century RSL change is constrained by long tide-gauge records from Pensacola, Florida (1925-present), Grand Isle, Louisiana (1947-present), and Galveston, Texas (1905-present) (see https://www.psmsl.org)(Fig. 4b).Pensacola is in a generally stable part of the coastline with negligible subsidence (Frederick et al., 2019): RSL rise for 1950-2010 was ~2.1 mm yr − 1 , but had accelerated to ~3.3 mm yr − 1 for 1993-2010.RSL rise at Grand Isle and Galveston differs from RSL rise at Pensacola due to subsidence, with RSL rise of ~9.3 mm yr − 1 at Grand Isle and ~ 6.8 mm yr − 1 at Galveston for 1950-2010.
In addition to secular RSL rise, MRD water-level fluctuations vary over a range of time scales and represent a variety of mechanisms (Hiatt et al., 2019).Diurnal tides are very small, with only ~330 mm in watersurface elevation change between daily mean high water (MHW) and mean low water (MLW) at Grand Isle (https://tidesandcurrents.noaa.gov/stationhome.html?id=8761724).However, tropical storm surges and river floods that coincide with cold front passage can raise water levels by up to 1 m (e.g., Roberts et al., 2015;Bilskie et al., 2016;Hiatt et al., 2019).Of importance to our study, Kolker et al. (2011) recognized annual anomalies in mean sea level (MSL) in the Pensacola, Grand Isle and Galveston tide-gauge time series, which are superimposed on the secular trend.Detrended annual anomalies in MSL for Pensacola, Grand Isle and Galveston for 1950-2010 are plotted in Figs.4c.The magnitude varies spatially along the coast, but anomalies covary between Grand Isle and Galveston with R = 0.84, and Grand Isle and Pensacola with R = 0.95, with P < 0.00001.The largest negative annual anomaly at Grand Isle was − 87 mm in 1963, whereas the largest positive annual anomaly was +96 mm in 1975, which sums to 183 mm of net RSL rise over that 12-year period, or ~ 15 mm yr 1 .By comparison, GMSL rise for the same period was a total of 7 mm, or 0.6 mm yr-1 (Dangendorf et al., 2017(Dangendorf et al., , 2019)).As discussed more fully below, the magnitude of positive annual anomalies in MSL can exceed rates of GMSL rise over decadal time scales by more than an order of magnitude, and the magnitude of negative annual anomalies in MSL can overwhelm rates of GMSL rise over decadal time scales to result in RSL fall.

MRD subsidence
Subsidence is intrinsic to large deltaic depocenters because of longterm crustal loading by sediment deposition.Long-term, time-averaged subsidence rates for the MRD vary from the north-to-south due to increases in thickness of the load from long-term deposition, but tend to be <− 0.5 to − 1 mm yr − 1 throughout the MRD (Fig. 5a) (Frederick et al., 2019).This long-term deep-seated component represents the accommodation that preserves deltaic deposits in the stratigraphic record, and is the sum of motion on growth faults, isostatic processes, and continued compaction of Pleistocene and older deposits.However, over late Holocene to modern time scales most subsidence has been generated by compaction of Holocene deltaic sediments (e.g., Törnqvist et al., 2008;Jankowski et al., 2017;Zoccarato et al., 2020).Compaction rates also increase with thickness of Holocene sediments and can be more than an order of magnitude higher than the long-term, time-averaged component associated with the deltaic depocenter.We note there are numerous deep-seated growth faults that have surface expression on the MRD delta plain (e.g., Gagliano et al., 2003;McCulloh and Heinrich, 2012;Shen et al., 2017;Frederick et al., 2019), and there is a clear spatial relationship between fault headwalls and land loss (Gagliano et al., 2003).However, rates and timing of fault motion are poorly constrained, and it is unclear whether fault movement causes land loss immediately afterwards, or whether fault slip at an earlier time created lower-elevation headwall topography that was preferentially submerged by sea-level rise at a later time (Roberts et al., 2008).
Subsidence rates for the period 1950-2010 can be estimated from differences between tide gauge records (e.g.National Research Council, 1987;Douglas, 2001;Zervas et al., 2013), where one of the gauges is assumed to be stable and undergoing minimal VLM.For example, subsidence at Grand Isle can be estimated by differencing Grand Isle RSL vs. Pensacola RSL, where VLM at Pensacola is less than ±1 mm yr − 1 (Fig. 4a): this approach yields an estimated mean subsidence rate of − 7.2 mm yr − 1 at Grand Isle for 1950-2010.Alternatively, Zervas et al. (2013) removed seasonal oceanographic effects (e.g. the inverse barometer effect and wind stresses), and differenced Grand Isle RSL vs. GMSL from Church and White (2011), and, to obtain an estimated mean subsidence rate of − 7.6 mm yr − 1 , whereas differencing Grand Isle RSL vs. GMSL from Dangendorf et al. (2019) yields an estimated mean subsidence rate of − 7.9 mm yr − 1 .GPS measurements provide a complementary regional perspective (e.g., Karegar et al., 2015;Byrnes et al., 2019;Byrnes, 2021;Wang et al., 2020;Hammond et al., 2021) but represent <20 yrs. of observations: in general, GPS vertical velocities are − 2 to − 3 mm yr − 1 to the north of ~30.2 • N within the MRD, and increase to − 5 to − 8 mm yr − 1 near the southern MRD shoreline at ~29 • N (Fig. 5b) (Byrnes et al., 2019;2021), which is generally consistent with subsidence rates estimated from differencing Grand Isle RSL and Pensacola RSL or GMSL.GPS vertical velocities are less than − 2 mm yr − 1 in locations along the Gulf of Mexico coast outside of the area of Holocene MRD deltaic sediments (Wang et al., 2020).Last, VLM has been estimated from GPS data for stations that are at, and/or surround, >2300 tide gauges worldwide (Hammond et al., 2021): vertical velocities at the Grand Isle GPS and 3 other nearby stations indicate VLM for the Grand Isle tide gauge is − 6.8 mm yr − 1 for 2006-2018.
An alternative subsidence dataset was published by Jankowski et al. (2017) based on the State of Louisiana's Coastwide Reference Monitoring System (LA-CRMS) network of almost 300 Rod Surface-Elevation Table and Marker Horizon (RSET) sites, where annual vertical accretion and surface-elevation change have been measured since the early 2000s, "shallow subsidence" is calculated as the difference between the two, and total subsidence is calculated as the sum of shallow subsidence plus vertical velocities derived from a north-south regression of GPS data (see Webb et al., 2013;Cahoon et al., 2019Cahoon et al., , 2020)).We discuss RSET vertical accretion rates below, but have concerns about the RSET method of calculating subsidence.For our analyses we use the subsidence model in Fig. 4b, which is derived from a regression of latitude vs. GPS vertical velocities that was initially published by Blum and Roberts (2012), and supplemented by new data from Kareger et al. (2015) and Byrnes et al. (2021).

MRD sediment supply and dispersal
Mississippi River suspended sediment loads have been measured at various locations since 1951 at or near Tarbert Landing, Mississippi (see Meade and Moody, 2010), which approximates the northern limits of the MRD.Numerous authors document the decline in loads that resulted from dams and other human activities in the Mississippi drainage basin, especially the closure of large dams in 1953 and 1955 on the Missouri River tributary some 2000 km upstream (e.g., Kesel et al., 1992;Blum and Roberts, 2009;Meade and Moody, 2010;Horowitz, 2010;Heimann et al., 2011aHeimann et al., , 2011b;;Kemp et al., 2016;Mize et al., 2018;Norton et al., 2019).Sediment loads at Tarbert Landing were ~ 460 Mt. yr − 1 (metric tons per year) for the three years of record before dam closure in but declined to <200 Mt. yr − 1 by 1980 and have continued to decline since (e.g., Mize et al., 2018;Norton et al., 2019).Recent estimates by Norton et al. (2019) further revised sediment loads downward to ~153 Mt. yr − 1 based on an evaluation and revision of collection methods and suspended sediment load calculations for Tarbert Landing.
Under natural conditions, sediment supplied to an active deltaic headland would be dispersed from the main distributary channel to the delta plain during flood events through a network of secondary channels.However, a continuous artificial levee system significantly limits riverine sediment dispersal regardless of the actual mass of sediment supplied to the delta region.By the mid 1800s, levees were generally continuous from just south of New Orleans to north of Baton Rouge on both sides of the Mississippi River, but a typical levee would have been only ~2.4 m high by the 1880s (USACE (United States Army Corps of Engineers), 2017).The modern levee system can be traced to large floods of the late 1800s, which prompted the United States government to construct a continuous line of defense against floods and increase recommended levee heights to ~6.7 m by 1914, to 8 m after the flood of 1927 and to 12 m after the flood of 1973 (Corthell, 1897;Smith and Winkley, 1996;Alexander et al., 2012; USACE (United States Army Corps of Engineers), 2017).Moreover, Bayou Lafourche, the primary distributary channel for the Lafourche deltaic headland when it was active prior to ~500 yrs.ago, remained partially connected to the Mississippi main stem but the remaining connection was sealed by a dam constructed in 1903 (Reuss, 2004).Mississippi River sediment loads have therefore been partially disconnected from the late Holocene deltaic headlands since the mid to late 1800s, and almost completely disconnected since the earliest 1900s.Recent studies show that much of the sediment that passes Tarbert Landing in the post-dam era is stored in the channel, between the levees, and does not reach the lowermost river south of New Orleans (Allison et al., 2012).
Longer-term sediment accumulation rates for the MRD can be estimated from radiocarbon dating of late Holocene flood-plain and deltaplain deposits.For example, Kesel (2008) obtained numerous radiocarbon ages from flood-plain deposits within the lower Mississippi Valley just to the south of Baton Rouge for the period ca.11-3.5 kyrs BP.Time-averaged rates of sediment accumulation were ~ 2.6 mm yr − over the entire period of record, but up to ~7 mm yr − 1 during initial development of the St. Bernard deltaic headland, and < 1 mm yr − 1 over the last 3.5 kyrs as the Lafourche and Plaquemines-Balize deltaic headlands prograded farther south.Similarly, Bomer et al. (2019) obtained radiocarbon ages from cores that sampled deltaic sediments to the south of New Orleans, on both the east and west side of the Mississippi River.These sites represent deposits of the St. Bernard, Lafourche, and Plaquemines-Balize deltaic headlands over the period ca.~3200-235 yrs.BP, and show that time-averaged sediment accumulation rates range from ~1 to 5 mm yr − 1 but are generally between 1 and 3 mm yr − 1 .A different perspective comes from crevasse splays of late Holocene age, or from measurements made after major flood events during the historic period, where rates of deposition can be several orders of magnitude higher than time-averaged values, and spatially variable.For example, Shen et al. (2015) shows that crevasse splay deposition rates can range from 10 to 40 mm yr − 1 , and can be sustained over centuries.Moreover, deposition during the 1973 flood ranged from 530 mm on natural levees to 11 mm in backswamp environments (Kesel et al., 1974), whereas deposition during the 2011 flood ranged from mm on natural levees to 3 mm in backswamps (Heitmuller et al., 2017).
Riverine sediment dispersal to the late Holocene deltaic headlands has been minimal since 1903, but the State of Louisiana's RSET stations, established in the early 2000s, show that organic-rich, low-density marsh sediment has accumulated on the southern Lafourche delta plain at rates of ~10-14 mm yr − 1 or more (Jankowski et al., 2017).Sanks et al. (2020) show the mass and volume of this marsh sediment is significant, and suggest the mineral component is delivered to the delta platform by tidal currents and storm-induced flooding, perhaps associated with the passage of cold fronts (e.g., Roberts et al., 2015).However, marsh sediments have bulk densities of 0.1-0.3g cm − 3 and organic fractions of 20-50% when deposited.By contrast, river-derived sediment measured by the RSET technique in the Plaquemines-Balize deltaic headland has significantly higher bulk densities (>0.6-0.9 g cm − 3 ) and smaller organic fractions (<15%) (Sanks et al., 2020), and late Holocene sediments of the Lafourche delta examined in Bomer et al. (2019) have bulk densities of >1 g cm − 3 after burial by ~250 mm or more of deposition, and organic fractions of <8%.

Prehistoric vs. historic period land loss
The delta cycle is autogenic and results in successive generations of deltaic headland construction that are lateral to each other, such that land gain and land loss occur simultaneously on active vs. abandoned deltaic headlands, respectively.Bentley et al. (2014) estimated that 40% of the MRD was growing through aggradation and progradation at any one time during the late Holocene, while 60% was undergoing submergence and degradation.They further argue that, over the middle to late Holocene, the length of deltaic coastline that was subject to land loss at any one time was always greater than the length of deltaic coastline that was experiencing land-area gain.The delta cycle is also timedependent to a considerable degree, and distinct MRD deltaic headlands are in different stages of degradation and submergence because they formed and were abandoned at different times (Gagliano and van   Beek, 1975).For example, pre-late Holocene deltaic headlands (the Maringouin and Teche deltas in Fig. 1) are known from subsurface data but have limited surface expression because they have degraded, subsided below sea level, and are now buried by late Holocene sediments (Boyd et al., 1989;Kulp et al., 2005).
There was net growth of the delta plain as a whole during the late Holocene, but land loss on individual deltaic headlands was ongoing prior to significant human intervention (see Boesch et al., 1994).In fact, as Gagliano et al. (1981) noted, a ~ 4000-year trend of net delta-plain growth was reversed by the late 19th century, with land loss emerging as a widely recognized existential crisis during the 20th century.Mid 1800s historic maps provide a qualitative perspective on the magnitude of pre-20th century land loss.Maps from Colton (1854), for example, show the St. Bernard deltaic headland, which formed ca.4000-1700 yrs.ago, already had a discontinuous delta plain that was partially submerged, and the delta front had been reworked into the Chandeleur barrier-island arc.Similarly, the Lafourche deltaic headland, which formed ca.1700-500 yrs.ago, retained primary distributary and interdistributary morphologies, but Terrebonne, Timbalier, and Barataria Bays had formed by conversion of distal interdistributary delta-plain environments to open water, and the delta front had been reworked into the Isle Derniers, Timbalier Island, and Grand Isle barrier-island arcs (Fig. 6).
Upper limits for the magnitude of pre-1932 land loss in the St. Bernard and Lafourche deltaic headlands can be approximated by comparing initial headland surface areas with land areas remaining in 1932 as measured in Couvillion et al. (2017) (Table 1).The St. Bernard headland is in the Pontchartrain and Breton Sound Coastal Hydrologic Basins, as defined by the State of Louisiana, which extend from Lake Pontchartrain to the Chandeleur barrier island arc (Figs. 2 and 6) (htt ps://lacoast.gov/new/About/Basins.aspx).These two basins have a surface area of ~14,000 km 2 , which we treat as a maximum value for the St. Bernard deltaic headland prior to abandonment ca.1700 yrs.BP.By 1932, land-surface area in these two coastal basins was ~4000 km 2 , hence up to ~10,000 km 2 had been converted to open water by that time, and an additional ~900 km 2 was lost between 1932 and 2016.Similarly, the Lafourche deltaic headland resides in the Barataria and Terrebonne Coastal Basins, with a combined surface area of ~14,000 km 2 , which we again treat as a maximum value for the Lafourche headland prior to its abandonment ca.~500 yrs.ago.By 1932, land area for the two basins was ~8300 km 2 , hence ~5700 km 2 had been converted to water by that time, and an additional ~2400 km 2 was lost between 1932 and 2016.From the above we estimate that as much as ~90% and ~ 70% of the land loss that has occurred in the St. Bernard and Lafourche headlands, respectively, occurred prior to 1932.

LiDAR digital topobathymetric modeling of historic-period land-area change
On a global scale, improved processing of digital-elevation data shows significantly more low-elevation coastal landscapes at-risk from global sea-level rise than previously thought (Kulp and Strauss, 2019).An integrated topographic and bathymetric dataset for the MRD region was published in, 2014 (CoNED, 2014: hereafter CoNED LiDAR): this dataset includes data collected from 2003 to 2013, with a spatial resolution of ~2 m, and a vertical uncertainty of ±20 cm.We use this integrated topobathymetric dataset to examine the land-loss issue in the southern MRD, which has experienced the highest 20th and early 21st century land-loss rates.Our focus area consists of ~14,500 km 2 , mostly within the Lafourche deltaic headland (see Fig. 7b), which was actively aggrading and prograding until ca.500 yrs.BP (see Chamberlain et al., 2018) and remained partially connected to the Mississippi River until 1903 (Reuss, 2004).  2 in Couvillion et al., 2017).Gray box represents the mean tidal range, from MHW to MLW.

Landscape sensitivity to small water-level changes
Elevation capital is a measure of marsh sustainability.Here, we use CoNED LiDAR to examine and quantify elevation capital for the MRD, and the sensitivity of land area to small changes in water levels of the scale produced by tides and/or annual anomalies in MSL.CoNED LiDAR is datumed to mean low water (MLW), equivalent to MSL-165 mm at the Grand Isle tide gauge in the southern Lafourche headland (https://tidesa ndcurrents.noaa.gov/datums.html?id=8761724).Hence, we created additional topobathymetric datasets that datum to MSL and to mean high water (MHW) at Grand Isle by adding 165 mm and 330 mm elevation, respectively, to each LiDAR data point.
We then used a binary classification of land (elevation ≥ MSL = 0) vs. water (elevation <0), which shows ~6350 km 2 is classified as land at MSL within our Lafourche deltaic headland study area.By comparison, Couvillion et al. (2017;hereafter Couv2017) classified ~6500 km 3 as land for 2014 within the same area, a difference of <3%.Our classification shows that >90% of the land-surface area in this part of the Lafourche deltaic headland is <0.5 m elevation: only distributarychannel alluvial ridges consistently exceed elevations of 0.5 m (Fig. 7a).Using the same classification approach, we calculated land areas of ~7580 km 2 with elevation ≥ MLW and ~ 4800 km 2 with elevation ≥ MHW (Fig. 7b), hence ~2780 km 2 of surface area is intertidal, and can be emergent or submerged during the normal tidal cycle.Last, we datumed CoNED LiDAR to Grand Isle MSL ± 80 mm, which simulates positive and negative annual anomalies in MSL, and find that ~1380 km 2 can be emergent or submerged with anomalies in MSL of this scale.
Fig. 7c shows that the largest measured departures in land area from Couv2017 for our Lafourche deltaic headland study area are ±350 km − 2 , which is a small value compared to land-area changes that can occur from the daily tidal cycle and/or anomalies in MSL.Metadata for Landsat images used in Couv2017 for land-area measurements was not available for pre-1985 images, but water levels at the Grand Isle tide gauge station were published for imagery used to measure land-area change from 1985 to 2016.Water levels departed from MSL by up to ±156 mm in any given year for images that cover our Lafourche deltaic headland focus area, and as shown in Fig. 8, there is an inverse relationship between water-level departures from MSL during image acquisition, and measured departures from the secular land-loss trend, with R = -0.81.In summary, we interpret the Couv2017 measurements to accurately represent land areas at the time measurements were made.However, the water surface was a non-stationary frame-of-reference during the period of measurements, and on a year-by-year basis the measurements themselves reflect the influence of changing water levels.

20th century land loss, subsidence, and sea-level rise
The century-scale secular land-loss trend has been attributed to varying degrees to ongoing subsidence and accelerated sea-level rise, coupled with reduced sediment dispersal to the delta plain.To further examine this issue, we hindcast land area at MSL, MHW and MLW for a series of time slices from 2010 to 1950 by progressively restoring landsurface elevation that has been lost because of subsidence and GMSL rise to the CoNED LiDAR dataset, while holding sediment input to zero.In other words, for each year in our hindcast we add the net elevation that has been lost from subsidence and GMSL rise to each LiDAR data point, then calculate land area with elevation ≥MSL.To estimate a range of land areas for each year in our simulation we also calculate land areas ≥MHW and ≥ MLW.
Our first hindcast uses steady subsidence rates derived from the latitude-dependent regression of GPS vertical velocity data in Fig. 5b, and values of GMSL rise of +0.7 mm yr − 1 for 1950-1970, +1.4 mm yr − 1 for 1971-1990, and + 2.5 mm yr − 1 for 1991-2010, as calculated from data in Dangendorf et al., 2019; see Fig. 9a).After restoring elevation lost to subsidence and GMSL rise for each annual time slice, individual pixels are again subject to a binary classification of land (elevations ≥0) vs. water (elevations <0) at MSL, MHW, and MLW.Results are then compared with measured changes in land area over time for this portion of the MRD from Couv2017.We also hypothesize that positive annual anomalies in MSL are especially significant to the inundation collapse of marsh environments, and increase land-loss rates because of their annual to multidecadal duration, whereas negative anomalies of the same magnitude and duration increase elevation capital and reduce land-loss rates.Hence, our second hindcast uses steady subsidence rates and rates of GMSL as above, but adds or subtracts anomalies in MSL from the Grand Isle tide gauge for each year in the simulation: we use conservative 10-yr anomalies in MSL following the view that it takes multiple years of flooding to cause marsh inundation collapse (e.g., Reed et al., 2020).
Results of our first hindcast are plotted in Fig. 9b, and show an accelerating rate of land loss from 1950 to 2010 that is well described by a 2nd order polynomial that mirrors the accelerating GMSL rise shown in Fig. 9a.Total simulated land loss is consistent with measurements: Couv2017 recorded land loss of 1860 km 2 between 1956 and 2010, which corresponds to a mean rate of ~34.5 km 2 yr − 1 , whereas our simulation produces 1833 km 2 for the same period, which corresponds to a mean rate of 33.9 km 2 yr − 1 .Our first simulation produced a total land loss of 1860 km 2 for the entire 1950-2010 time period, which corresponds to ~31 km 2 yr − 1 , and the Couv2017 measurements for all years reside between our simulated land areas at MHW and MLW.However, the 2nd order polynomial fit overestimates total land area for the 1980s and 1990s, and overestimates land-loss rates after 1995.Results from our second hindcast are plotted in Fig. 9c, and are again consistent with Couv2017, but produce a closer correspondence between simulated values and measurement for most individual years.Incorporation of anomalies in MSL into the hindcast modestly increased land loss rates for 1950-2010 as a whole to 31.8 km 2 yr − 1 , and all measurements from Couv2017 again reside within the bounds set by simulated land areas at MHW and MLW for that year.Inclusion of anomalies in MSL also results in a 4th order polynomial fit that features acceleration of land-loss rates in the late 1960s and 1970s, then deceleration in the late 1980s and 1990s.

How much riverine sediment was needed to mitigate 20th century land loss?
In a 2-dimensional north-south profile through our Lafourche deltaic headland focus area, the sediment thickness that would have been needed to mitigate land loss is the inverse of accommodation generated by subsidence and GMSL rise.Hence, we use CoNED LiDAR to model the riverine sediment input that would have been required to fill accommodation that was generated and mitigate the observed trend in land loss.For 1950-2010, a period over which GMSL rise accelerated from ~0.7 mm yr − 1 to ~2.5-3 mm yr − 1 (e.g., Dangendorf et al., 2019), this would result in accommodation of ~3-4 mm yr − 1 to the north, where subsidence rates are low and accommodation is dominated by GMSL rise, vs. 9-10 mm yr − 1 for the southernmost deltaic shoreline where subsidence rates are − 6 to − 8 mm yr − 1 and subsidence dominates accommodation.For comparison, these values are equivalent to spatially averaged means of ~5.5 mm yr − 1 of net deposition with zero compaction in 1950 when GMSL rise was ~1.5 mm yr − 1 , and ~ 7 mm yr − 1 of net deposition with zero compaction in 2010 when GMSL rise was ~3 mm yr − 1 .
Rates of deposition that would have been required to mitigate land loss are more than twice as high as time-averaged sediment accumulation rates based on radiocarbon dating of the prehistoric stratigraphic record, a deficit consistent with the century-scale secular trend of land loss.But they also raise questions about the significance of the organicrich, low-density sediment that has been documented using the LA-CRMS RSET sites.This component of deltaic sedimentation has likely always been present in the Lafourche deltaic headland focus area and elsewhere, not just during the ~15 yrs.over which RSET measurements have been collected.However, it has likely been the only significant component of deltaic sedimentation after the connection between the Mississippi River and Bayou Lafourche was severed in 1903.Vertical accretion of this low-density marsh sediment has proceeded at rates of 10-14 mm yr − 1 in the southern Lafourche deltaic headland during the period of RSET data collection, which exceeds the generation of accommodation from subsidence and sea-level rise by a factor of 2. However, Couv2017 shows the Lafourche deltaic headland has sustained the highest rates of land loss in the MRD over this same time period.
These results should be viewed within the context of the role of compaction in sustaining surface elevation.As noted above, rapid compaction of shallow mud-and organic-rich deltaic sediment is the primary driver for subsidence over time scales of decades to millennia (e.g., Törnqvist et al., 2008), and therefore plays a strong role in the net vertical accretion of the land surface over those time scales.Assume, for example, that compaction will, over decades to centuries, increase the bulk density of newly deposited sediment to a value that approaches typical shallow subsurface bulk densities of 1.5 g cm − 3 .Conceptually, marsh sediment with an initial bulk density of ~0.3 g cm − 3 will eventually compact to ~20% of its original thickness, whereas riverine sediment with an initial bulk density of ~1 g cm − 3 will eventually compact to 67% of its original thickness.In terms of sustaining landsurface elevation over decadal to century time scales, 5.5 mm yr − 1 of spatially averaged sediment accretion without compaction in 1950 corresponds to ~27.5 mm yr − 1 of organic-rich marsh sediment with an initial bulk density of 0.3 g cm − 3 vs.~8.2mm yr − 1 of riverine sediment deposition with an initial bulk density of 1 g cm − 3 .

General considerations
The MRD is a vast wetland landscape with minimal elevation capital, and is sensitive to the balance between deposition of sediment and loss of elevation due to relative sea-level rise.>90% of the land area in our Lafourche deltaic headland study area is <50 cm in elevation relative to MSL, and almost 2800 km 2 can be submerged or emergent from the daily tidal cycle.Moreover, almost 1400 km 2 can be submerged or emergent due to annual anomalies in MSL of ±80 mm.We view annual anomalies in MSL to be especially important to the land loss issue because they raise or lower average water levels for multiple years at a time: as argued in Reed et al. (2020), MRD marshes are resilient with respect to shortterm flooding from tidal cycles, but inundation collapse and conversion of marsh to open water can result from flooding of marshes for as little as two years.Annual anomalies in MSL therefore provide an actual mechanism for acceleration and deceleration of land loss on multiannual to multidecadal time scales.
The empirical trend of land loss in the MRD for 1932-2016 is well documented in Couv2017.Our first hindcast using CoNED LiDAR for 1950-2010 focused on land-area change due to steady subsidence and accelerating GMSL rise, and produced values for total land loss and timeaveraged rates that closely correspond to values from Couv2017.However, this hindcast shows a secular trend of continually accelerating land-loss rates over time that mirrors the acceleration of GMSL rise (Fig. 9b).Our second hindcast incorporates GMSL plus 10-yr anomalies in MSL from Grand Isle, and also produces values for total land loss and time-averaged rates that closely resemble values from Couv2017, but features acceleration of land loss in the late 1960s and 1970s, and deceleration in the mid-late 1990s Fig. 9c).As noted above, there has long been discussion of the relative influence of the natural delta cycle vs.human activities on MRD land loss (e.g., Penland et al., 1990;Boesch et al., 1994;Kolker et al., 2011).Our results indicate that the general trajectory and magnitude of 20th and early 21st century land loss reflects accelerating GMSL rise with reduced riverine sediment input on a continually subsiding delta plain, and consideration of other factors is not necessary to account for observed land loss.Inclusion of 10-yr anomalies in MSL from Grand Isle in our second hindcast only changes the total land loss over the period 1950-2010 by <1 km 2 yr − 1 , with the effects of positive vs. negative anomalies essentially canceling each other over the 1950-2010 period of our study.
A number of studies have attributed acceleration and deceleration of land-loss to non-linear VLM that reflects subsurface fluid extraction (e. g., Morton et al., 2003Morton et al., , 2005Morton et al., , 2006;;Mallman and Zoback, 2007;Chan and Zoback, 2007;Kolker et al., 2011;Chang et al., 2014;Day et al., 2020).We propose a different interpretation.First, land-loss departures from a linear secular trend, as measured in Couv2017, are only <20% and < 50% of the land-area change that can result from the daily tidal cycle or annual anomalies in MSL, respectively (Fig. 7).Second, there is a clear inverse relationship between measured water-level departures from MSL at Grand Isle for the time of acquisition of imagery used to measure land-area change, and the departures in land loss from the secular trend (Fig. 8).Third, all Couv2017 measurements for 1950-2010 reside between our hindcasted land areas at MLW and MHW if we use steady subsidence and GMSL rise alone, or when we include 10-yr anomalies in MSL (Fig. 9).We conclude that observations of acceleration and deceleration of land-loss rates should be taken at face value, but they represent measurements that were made when water levels departed from MSL due to the tidal cycle and/or annual anomalies in MSL.
Our topobathymetric hindcasts held riverine sediment input to zero.However, land loss from 1950 to 2010 could have been mitigated if riverine sediments had been dispersed to the delta plain at rates equal to the accommodation produced by relative sea-level rise.From a massbalance perspective, these values are equivalent to a spatially averaged rate of ~5.5 mm yr − 1 in 1950, accelerating to 7 mm yr − 1 by 2010, and are generally higher than, but within the range of, time-averaged rates of deposition for riverine sediments that have been measured from the late Holocene to modern stratigraphic record.However, interestingly, these values are <50% of the vertical accretion rates for organic-rich, low-density marsh sediment that have been measured at the State of Louisiana's RSET sites within our Lafourche deltaic headland over the last 2 decades.Such organic-rich, low-density marsh sediment accumulation has likely been the norm since the connection between Bayou Lafourche and the Mississippi River was severed in 1903 (Reuss, 2004), yet the Lafourche headland has sustained the highest 20th century MRD land loss rates.As noted in Keogh et al. (2021), if sediment is introduced in the delta region to sustain surface elevation, the elevation gained from new sediment deposition must outpace the elevation lost to compaction.We question whether this will ever be the case for the organic-rich low-density marsh sediment measured at RSET sites, and suggest that, over multidecadal and longer timescales, it compacts rapidly and, on a year-to-year basis, plays an insignificant role in maintaining or increasing surface elevation.Modeling of compaction by Zoccarato et al. (2020) and Keogh et al. (2021) supports the view that most of the marsh sediment thicknesses recorded at RSET sites will be lost to compaction over multidecadal and longer time scales.
In summary, we view MRD land-loss through the lens of the delta cycle, where ongoing subsidence, reduced riverine sediment input, and acceleration of GMSL rise promotes rapid degradation and submergence, shoreline transgression, and landward migration of brackish and saline marsh environments.We also view the land-loss issue through the lens of the Holocene stratigraphic record.Törnqvist et al. (2020) show that marsh deposits in subsurface cores are commonly overlain by lagoonal muds, indicating they are part of an overall shoreline transgression, and interpret the presence of marsh deposits to indicate that a tipping point had been reached and submergence was inevitable.It is worth noting that marsh drowning occurred within a few centuries of marsh formation when time-averaged rates of RSL rise, measured from the stratigraphic record as well, exceeded ~3 mm yr − 1 , but when timeaveraged RSL rise exceeded 6 to 9 mm yr − 1 , marsh environments were converted to open water in ~50 years.

Causes of anomalously high or low land-loss rates
20th and early 21st century land-loss trends reflect an unfortunate convergence of ongoing subsidence, accelerating GMSL rise, and limited riverine sediment dispersal to the delta plain due to levee construction, and we show that land-loss departures from a secular linear trend reflect measurements taken when water levels departed from MSL due to the tidal cycle, annual anomalies in MSL, or both.We view positive anomalies in MSL to be especially important for acceleration of land-loss rates because they can increase water depths on marsh platforms for multiple years at a time, and can therefore trigger inundation collapse (e.g., Reed et al., 2020).By contrast, negative anomalies in MSL can decrease water levels on marsh platforms for multiple years at a time, and can result in deceleration of land loss from inundation collapse.Annual anomalies in MSL are mapped from the global network of tide gauges, the effects of anomalies on estuarine and coastal processes are becoming widely appreciated, and anomalies have been linked to atmospheric and oceanographic drivers that operate 100 s to 1000s of kilometers distant (e.g., Nerem et al., 1999;Morris, 2000;Kolker et al., 2009;Sweet and Zervas, 2011;Merrifield et al., 2012;Thompson and Mitchum, 2014;Theuerkauf et al., 2014;Hamlington et al., 2015Hamlington et al., , 2020;;Barnard et al., 2017;Sweet et al., 2018;Han et al., 2017Han et al., , 2019;;Volkov et al., 2019;Woodworth et al., 2019;Chafik et al., 2019).Examples of annual anomalies in MSL for North America are illustrated in Fig. 10.Kolker et al. (2011) first recognized regional coherence in annual anomalies in MSL at the Pensacola, Grand Isle, and Galveston tide gauges, and attributed them to reflect annual wind and sea-surface temperature anomalies.Kolker et al. (2011) then calculated a linear subsidence rate of − 7.6 mm yr − 1 for Grand Isle from 1947 to 2006 by differencing the Grand Isle vs. Pensacola tide gauge time series, which is within the range of results obtained by other workers (e.g., Zervas et al., 2013; our calculations reported above) over the years.However, they then assumed that residuals in the difference between Grand Isle and Pensacola reflected changes in VLM from changes in rates of subsurface fluid withdrawals.They calculated subsidence rates of ~ − 3.2 mm yr − 1 for the first decade of the record, ~ − 12.6 mm yr − 1 for 1959-1974, and ~ − 1 mm yr − 1 for 1992-2006.We disagree with the interpretation of these residuals as VLM, as follows: 1. Long tide-gauge records from around the GoM show the magnitude of annual anomalies in MSL varies by a factor of 2, and differences in the magnitude of anomalies between stations is the norm.Anomalies are larger to the north and west, and smaller along the Florida peninsula to the east and south (Fig. 11a).It is not clear why Pensacola would represent a norm for comparing the magnitude of anomalies at other stations.2. The interpreted period of maximum subsidence at Grand Isle in Kolker et al. (2011) does correspond to high rates of hydrocarbon production (see Day et al., 2020), but it also corresponds with a multi-decadal transition from generally negative (1950s and 1960s) to generally positive (1970s to 1990s) annual anomalies in MSL at Grand Isle, as well as at tide-gauges across the GoM and along the SE Atlantic coast (Fig. 10).3.If interpretations of non-linear VLM are correct, the secular subsidence rate for 1950-2010 should resemble the low rates estimated for years prior to and after the inferred effects of hydrocarbon production, before 1968 and after 1995, which would be <− 3.2 mm yr − 1 at Grand Isle.We regard such low secular subsidence rates to be unlikely.4. Campaign-style GPS data collected from 2003 to 2018 (e.g., Byrnes et al., 2019;2021) show subsidence rates of − 6.0 to − 8.0 mm yr − 1 at the Grand Isle GPS station and other sites close to the deltaic shoreline (see also Hammond et al., 2021).These rates are consistent with the 1950-2010 secular trend from differencing the Grand Isle and Pensacola tide gauge records, but inconsistent with the low rates calculated by Kolker et al. (2011) for the immediately preceding 1993-2006 period.
We therefore interpret positive and negative anomalies in MSL at Grand Isle to be larger in magnitude than those at Pensacola and elsewhere in the GoM, and generally unrelated to VLM.
This interpretation begs the question of what causes the anomalies, and why are they larger at Grand Isle than elsewhere?We do not seek to delineate all factors that contribute to anomalies in MSL at Grand Isle, but we briefly discuss drivers that have been featured in the literature.First, at the largest scale, Thompson and Mitchum (2014) and Volkov et al. (2019) show that annual anomalies in MSL from the GoM covary with those of the SE US Atlantic coast to the south of Cape Hatteras, North Carolina, and that Atlantic water advected into and out of the GoM likely plays a significant role in generating positive and negative anomalies in the GoM.Volkov et al. (2019) also show that the leading mode (empirical orthogonal function 1) in Atlantic subtropical gyrescale variability in sea surface height is responsible for most of the interannual-to decadal-scale sea surface height changes along the SE US Atlantic and the Gulf of Mexico.They calculate a correlation coefficient of 0.69, significant at the 95% level, between changes in the leading mode over time and the detrended Grand Isle tide gauge record, and discuss how the leading mode is related to large-scale atmospheric circulation and the Atlantic Meridional Overturning Circulation (AMOC).
On the more local scale, early work by Meade and Emery (1971) shows that variations in river discharge to the coastal ocean can also drive annual anomalies in MSL.Discharge from river mouths produces a lens of freshwater that floats on the more dense saltier water of the coastal ocean (Durand et al., 2019;Hamlington et al., 2020).Freshwater discharge plumes are then subject to Coriolis deflection (e.g., Piecuch et al., 2018;Fofonova et al., 2021), which, for the northern GoM, results in transport of freshwater alongshore and to the west (Fig. 11c and d).
Freshwater discharge is especially pertinent to annual anomalies in MSL at Grand Isle because of its location between the Mississippi and Atchafalaya River mouths, the two discharge points for the Mississippi River system, which accounts for ~90% of freshwater discharge to the northern GoM.Meade and Emery (1971; see also Piecuch et al., 2018), concluded that annual anomalies in MSL of ± ~0.1 mm can be expected for each km 3 in anomalous annual discharge to the Gulf of Mexico, "if spread evenly over the surface".However, neither the Meade and Emery (1971) or Piecuch et al. (2018) studies included data from Grand Isle, which has the largest annual anomalies in MSL for the Gulf of Mexico.Regardless, mean annual discharge for the Mississippi River system (the combined Mississippi and Atchafalaya Rivers) is ~676 ± 184 km 3 yr − 1 for 1950-2010 (from data in Piecuch et al., 2018), and annual departures approach ±400 km 3 yr − 1 (see Fig. 12c).From the basin-scale relationships in Meade and Emery (1971) and Piecuch et al. (2018), annual discharge anomalies of this scale can produce annual basin-scale MSL anomalies of ± ~40 mm, which is close to the maximum difference between annual anomalies in MSL at Pensacola and Grand Isle.We suggest that basin-scale values such as these represent minimum values for freshwater discharge influences on the Grand Isle tide gauge.
The above drivers change water flux into the Gulf of Mexico, and are related to larger-scale atmospheric drivers.We therefore followed up on the Kolker et al. (2011) argument for atmospheric influence on anomalies, and used the European Center for Medium-Range Weather Forecasting (ECMWF) Reanalysis version 5 (ERA5; Hersbach et al., 2020) to obtain annual mean sea-level pressure, 10-m wind, and precipitation anomalies in the Gulf of Mexico watershed for years that are within the upper and lower quartiles of annual anomalies in MSL at Grand Isle (Fig. 12a).Similar to Kolker et al. (2011), we find the upper quartile corresponds to anomalously high pressure over the subtropical North Atlantic and low pressure over the northcentral US (Fig. 13a).This pattern enhances moisture flux into the US continental interior and Gulf of Mexico watershed, which leads to positive anomalies in precipitation and Mississippi River discharge (Figs.12b, c and 13b).Strong SSE winds on the western and southern side of the MRD, and associated Ekman transport in the shallow shelf waters, then pile up water along the Lafourche deltaic headland coast to produce higher coastal sea level and higher rates of land loss.By contrast, the lower quartile corresponds to anomalously high surface pressure and anticyclonic flow over the western Gulf of Mexico (Fig. 13c), with reduced moisture flux into the continental interior, and negative anomalies in precipitation and Mississippi River discharge (Figs.12b, c and 13d).Weak northerly winds along the northern Gulf coast then favor Ekman transport offshore, which lowers coastal sea level and reduces land loss.Our proposed interpretation therefore includes advection of Atlantic water into and out of the GoM, amplification of positive and negative anomalies in Grand Isle MSL due to anomalies in freshwater discharge from the Mississippi River system, and further amplification of positive and negative anomalies at Grand Isle and the Lafourche deltaic headland due to strong SSE vs. weak northerly winds.
The multidecadal component of land loss is critical to address as well.As shown in Figs. 12, 10-yr anomalies in Grand Isle MSL, precipitation over the Gulf of Mexico watershed, and discharge from the Mississippi River system were below the secular trend prior to the midlate 1960s, above the secular trend from the late 1960s to the mid-late 1990s, and then below the secular trend until 2010.10-year anomalies in Grand Isle MSL are positively correlated with 10-yr anomalies in Mississippi River discharge, with R = 0.73, and P < α = 0.01.Moreover, 10-yr anomalies in MSL are sufficient in magnitude to overwhelm GMSL rise and the secular rates of subsidence to change the magnitude and direction of RSL change, as well as cause acceleration and deceleration in MRD land loss rates over multidecadal time scales.For example, in Fig. 9a we add 10-yr anomalies in Grand Isle MSL to GMSL for use in our hindcast shown in Fig. 9c.From 1951 to 1965, a period characterized by negative 10-yr anomalies in MSL, 10-yr RSL fell at − 1.5 mm yr − 1 , and remained 20-40 mm below GMSL until 1966, which is consistent with low rates of land loss for that time period.By contrast, 1966By contrast, -1973 represents the transition from negative to positive 10-yr anomalies in MSL: 10-yr RSL rose at 12.6 mm yr − 1 until 1973, then slowed to 0.3 mm yr − 1 but remained 20-40 mm above GMSL until 1994, which is consistent with acceleration of land loss rates from marsh inundation collapse.Beginning in 1994, when negative anomalies in MSL again became the norm, 10-yr RSL fell at − 3.1 mm yr − 1 until 2002 but then rose 2.6 mm yr − 1 from 2002 to 2009 and remained ~20-30 mm below GMSL through 2009, which is consistent with observations of deceleration of land loss and lower land-loss rates.
Several oceanographic and climate cycles have been associated with anomalies in MSL, but the Atlantic Multidecadal Oscillation (AMO) provides the greatest explanatory power at Grand Isle for 1950-2010.The AMO is described empirically by slowly fluctuating sea-surface temperature (SST) anomalies in the North Atlantic, and has positive (warm) and negative (cool) phases that average 30 ± 5 years in duration (Schlesinger and Ramankutty, 1994;Kerr, 2000;Trenberth and Shea, 2006).Moreover, it has been suggested that the AMO is coupled to, and covaries with, the Atlantic Meridional Overturning Circulation (AMOC) (Trenary and DelSole, 2016;Kim et al., 2021;although see Clement et al., 2015see Clement et al., , 2016;;Zhang et al., 2016), which has been associated with transport of Atlantic water into and out of the GoM, as discussed above (e.g.Thompson and Mitchum, 2014;Volkov et al., 2019;Han et al., 2019).Previous work shows the AMO correlates to a range of atmospheric effects (e.g., Enfield et al., 2001;Knight et al., 2006), including droughts and wet periods over the US Great Plains (Hu et al., 2011;Nigam et al., 2011;Hu and Feng, 2012).AMO warm phases of the 1950s to mid-1960s, and then after 1995, are associated with droughts in the GoM watershed, negative anomalies in Mississippi River discharge, negative anomalies in MSL at Grand Isle, and lower measured land loss.By contrast, the AMO cool phase of the mid-1960s to the mid-1990s is associated with positive anomalies in precipitation in the GoM watershed, positive anomalies in Mississippi River discharge, positive anomalies in Grand Isle MSL, and higher measured land loss.10-year anomalies in Mississippi River discharge and Grand Isle MSL are inversely correlated to the AMO index, with R = − 0.56 and R = − 0.86, respectively, with P < α = 0.01.
There are different views about whether the AMO reflects a natural oscillation in the atmospheric and oceanographic system, or whether it reflects natural and anthropogenic forcing (see Mann et al., 2014Mann et al., , 2020Mann et al., , 2021;;Müller-Plath, 2020).Regardless, the empirical AMO index describes slow variations in the atmosphere-ocean circulation for 1950-2010.Our conceptual model is based on relationships between advection of Atlantic water into and out of the GoM, and the AMO's impact on precipitation in the GoM watershed, Mississippi River discharge, anomalies in MSL at Grand Isle, and MRD land-loss rates (Fig. 12).These relationships provide insight into variations on time scales between the year-to-year variability expected from annual anomalies in MSL alone and the century-scale secular trend, and they highlight the importance of this slow oscillation in modulating the effects of relative sea-level rise in the northern GoM and land-loss rates in the MRD.

Conclusions
The Mississippi River delta suffered large-scale land loss during the 20th century and is representative of a global-scale phenomenon where low-elevation deltaic coasts are increasingly at risk from sea-level rise and disruption of the sediment-dispersal system.The late Holocene MRD was constructed when GMSL rise was <1 mm yr − 1 : deltaic landscape evolution followed a cyclical pattern where an active deltaic headland aggraded and prograded seaward, then the active river channel avulsed and relocated to construct a new headland laterally along the coast.The abandoned delta front was then eroded by waves and currents and the abandoned delta plain subsided, degraded, and was ultimately submerged.Maps from the mid-1800s show the St. Bernard deltaic headland, which formed ca.4000-1700 yrs.BP, had significantly degraded by that time, and no longer retained primary morphological characteristics, whereas the Lafourche deltaic headland to the west and south, which formed ca.1700-500 yrs.BP, retained primary morphology but the land-loss trend was already established in low-elevation interdistributary areas (see also Twilley et al., 2016).The modern Plaquemines-Balize "birdsfoot" deltaic headland was well developed by this time.
Disruption of sediment-dispersal from engineering activities, especially levee construction, increased through the 1800s to the point where Corthell (1897) recognized that a continuous artificial levee system along the Mississippi River would limit sediment dispersal to the broader delta plain, which would in turn eventually result in subsidence of "delta lands below the level of the sea".Collectively, river engineering had significantly disrupted sediment dispersal by the early 1900s, which likely accelerated land-loss rates as Corthell (1897) had predicted.However, unlike the effects of levee construction, the 20th century acceleration in GMSL rise was not foreseen, and by the end of the 20th century rates of GMSL rise were 3-4 times as fast as during the late Holocene when the St. Bernard and Lafourche deltaic headlands were constructed.Land loss in coastal Louisiana is an existential crisis, and continued land loss is predicted for all environmental scenarios by 2060 (Reed et al., 2020), with near complete submergence of the MRD by 2100 (Blum andRoberts, 2009, 2012; see also Kulp and Strauss, 2019).
Multidecadal periods of land-loss acceleration and deceleration are superimposed on the century-scale trend (Couvillion et al., 2017).We investigated the sensitivity of the MRD to changes in water levels, as well as the trajectory of land loss for 1950-2010, and find the following: • The MRD has inherently limited elevation capital due to the low tidal range, and >90% of the surface area in the Lafourche deltaic headland is <0.5 m elevation.Moreover, ~2800 km 2 of surface area in our study area can be submerged or emergent depending on the tidal cycle, and ~ 1400 km 2 can be submerged or emergent from annual anomalies in MSL.The widely discussed departures from the secular trend of land loss in our study area, interpreted to represent accelerations and decelerations in land loss rates, are small by comparison, only ±300 km 2 .We conclude these departures reflect measurements that were made when water levels departed from mean sea level due to position within the daily tidal cycle and/or anomalies in MSL.• Total MRD land loss and time-averaged land-loss rates for 1950-2010 are adequately explained by ongoing subsidence and GMSL rise, in the absence of significant sediment input, but these variables alone do not produce observed periods of land loss acceleration and deceleration.If we also consider 10-yr anomalies in MSL, our hindcasted land areas again closely correspond to measured land loss and time-averaged rates, but show acceleration during the late 1960s and 1970, and deceleration during the late 1980s and 1990s.We conclude that annual and 10-yr anomalies in MSL are of sufficient magnitude to overwhelm GMSL and change the magnitude and direction of RSL, as well as modulate land loss rates over annual to multidecadal timescales.But the net effects of positive and negative anomalies over 1950-2010 are very small because they mostly cancel each other out over that time period, and therefore do not significantly change the secular trend in land loss.• Our estimates for how much riverine sediment deposition would have been sufficient to maintain a steady-state land area for the Lafourche headland in 1950 are ~5-7 mm yr − 1 , which are within the range of late Holocene and modern deposition rates.But they are also <50% of typical vertical accretion rates for the organic-rich lowdensity marsh sediment that have been measured for the Lafourche headland, which has experienced the highest land loss rates in the MRD.We conclude that this apparent contradiction reflects the ineffectiveness of low-density organic-rich marsh sediment to sustain surface elevation because of high and rapid compaction rates.
In summary, the 20th and early 21st century secular trend of MRD land loss reflects an unfortunate convergence of a subsiding delta plain, acceleration of GMSL rise, and greatly reduced sediment dispersal from levee construction.Multidecadal changes in land-loss rates are superimposed on the secular trend, and are caused by naturally occurring, regionally coherent annual to multidecadal anomalies in MSL (Fig. 14).Anomalies in MSL, in turn, reflect an ensemble of processes that begins with large-scale transport of water into and out of the Gulf of Mexico, and temperature anomalies in the North Atlantic Ocean that are associated with the Atlantic Multidecadal Oscillation, which then drives anomalies in precipitation in the Gulf of Mexico watershed, Mississippi River system discharge, and coastal wind stresses.The annual to multidecadal anomalies in MSL described here are well documented on a global scale and will modulate the secular trend of GMSL rise and its impacts on delta plains elsewhere as they respond to ongoing subsidence and anthropogenic disruption of sediment dispersal.

Fig. 1 .
Fig. 1.(a) General location of the Mississippi River Delta study area, illustrating the Mississippi River contributing drainage basin.Red box indicates location of (b) and (c).(b) Generalized geologic map illustrating distribution of late Pleistocene and Holocene environments of deposition.(c) 2014 Landsat image with polygons that illustrate the general location and extent of major deltaic headlands that comprise the MRD.Modified from Blum and Roberts (2012) and references therein.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Mississippi River Delta land loss from 1932 to 2016.(a) Landsat 8 image showing the generalized extent of Late Holocene deltaic headlands (after Blum and Roberts, 2012 and based on Frazier, 1967).Areas of 1932-2016 land loss are shown in red and areas of land gain are shown in yellow.New Orleans is designated by NOLA, whereas Grand Isle, Louisiana is designated by GI.(b) Land-area measurements from 1932 to 2016 (black dots) for the Mississippi River delta as a whole, including the Louisiana Chenier Plain to the west, with a spline best fit model (red line) and 95% confidence bands (thin dashed black lines).(c) Modelled trend of land-loss rates based on a spline best fit model (red line) and confidence bands (thin dashed black lines).All land-loss data are from Couvillion et al. (2017).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 4a plots GMSL from Dangendorf et al. (2019) for our 1950-2010 period of interest, illustrating the rates we use in our analyses.Relative sea-level (RSL) change for the northern GoM over multidecadal time scales is a function of GMSL change and other factors,

Fig. 4 .
Fig. 4. (a) Rates of GMSL rise used for simulations in Fig. 9b and c (from Dangendorf et al., 2019).(b) Relative sea-level change at the Grand Isle (Louisiana), Galveston (Texas) and Pensacola (Florida) tide gauges, with best-fit linear regression models (dashed lines).Each time series has been normalized to a value of 0 for 1950 to enable direct comparison.Note that Grand Isle is missing annual MSL data from 1976 and 1977.Data were obtained from the UK Permanent Service for Mean Sea Level at https://www.psmsl.org.(c) Detrended annual anomalies in MSL from the Grand Isle, Galveston, and Pensacola tide-gauge records.

Fig. 5 .
Fig. 5. (a) Time-averaged long-term vertical land motion characteristic of the MRD deltaic depocenter for a north-south profile through the Lafourche deltaic headland (from Frederick et al., 2019).General latitudes of Baton Rouge, New Orleans, the Empire-Golden Meadows fault zone (EGMFZ), and Grand Isle as shown.(b) Early 21st century vertical velocities and uncertainties from campaign-style GPS stations within the St. Bernard and Lafourche deltaic headlands plotted against latitude (from Byrnes et al., 2019; 2021).Depth of penetration for mounting rods or structures as shown.

Fig. 6 .
Fig.6.Early historic map byColton (1854), which illustrates the different stages of degradation of the St. Bernard and Lafourche deltaic headlands by the mid 18th century.Dashed red lines define the limits of the Ponchartrain and Breton Sound coastal basins, as defined by the State of Louisiana (https://lacoast.gov/new/About/Basins.aspx),which we use as a proxy for the maximum surface area of the St. Bernard deltaic headland, and the Barataria and Terrebonne coastal basins, which we use as a proxy for the maximum surface area of the Lafourche deltaic headland.The seaward limits of these coastal basins generally approximate the maximum deltaic shorelines estimated from the extent of delta-front sand byFrazier (1967).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. (a) Frequency distribution of elevation in the southern MRD from 2014 CoNED LiDAR topobathymetric data, measured as the number of LiDAR datapoints within the designated elevation range.(b) 2014 Landsat 8 image of our Lafourche deltaic headland focus area, with simulated land-area change at MHW and MLW as shown.(c) Summary of land-area change for the Lafourche deltaic headland at MHW and MLW, from annual anomalies in MSL of ±80 mm, and from annual anomalies in MSL ± MHW and MLW.The gray box shows the extent of measured departures in land loss from the linear secular trend for our Lafourche deltaic headland focus area, which have been used to infer acceleration and deceleration of land loss over time (from Couvillion et al., 2017).

Fig. 8 .
Fig. 8. Inverse relationship between measurements of land loss and water levels at the time of image collection (from Couvillion et al., 2017), plotted as a time series (a) and as a scatter plot (b).Land loss departures represent departures from the secular trend for the Lafourche deltaic headland study area in Fig. 7b.Water level departures represent the average water levels at the time of collection of Landsat images, relative to Grand Isle MSL (from Table2inCouvillion et al., 2017).Gray box represents the mean tidal range, from MHW to MLW.

Fig. 9 .
Fig. 9. (a) Plot of GMSL rise from Dangendorf et al. (2019), and GMSL ±10-yr anomalies in Grand Isle MSL.(b) Hindcasted land area at MSL for 1950-2010 in response to subsidence from Fig. 5b and GMSL rise from Fig. 9a, with sediment input = 0, compared to land area measurements from Couvillion et al. (2017).(c) Simulated land area at MSL through time in response to subsidence from Fig. 5b and GMSL rise ± Grand Isle sea-level anomalies from Fig. 9a, with sediment input = 0, compared to measurements from Couvillion et al. (2017).Dashed blue lines represent the best fits for our simulated values, whereas vertical lines with endcaps represent the range of land area values between MLW and MHW.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Regional-scale annual anomalies in MSL for the years indicated from tide-gauge data in North America.(a to d) Negative annual MSL anomalies in the Gulf of Mexico.(e to h) Positive annual anomalies in MSL in the Gulf of Mexico.Each map represents departures for that tide-gauge station from the mean for the period 1960-1990.Maps obtained from the UK Permanent Service for Mean Sea Level (PSMSL) at https://www.psmsl.org.

Fig. 11 .
Fig. 11.(a) Box and whisker plot of annual anomalies in MSL from Gulf of Mexico tide gauges, illustrating regional variability in the magnitude of anomalies, with the highest values at Grand Isle, Louisiana, and the lowest values at Key West, Florida.The box defines the middle two quartiles of the distribution, whereas the horizontal line defines the median, and data points outside the whiskers are statistical outliers.(b) Conceptual model for freshwater discharge and generation of a freshwater lens at the Grand Isle tide gauge in streamwise cross-section, where freshwater floats on saline ocean waters, and forms a bulge, or lens, at the top of the water column.(c) Conceptual model for freshwater discharge and generation of a freshwater lens in map view, where freshwater is steered to the west by Coriolis deflection, which provides an alongshore component.In (b), h w = water column height in the river channel, and h lens = thickness of the freshwater lens in the nearshore GoM.In (c), width of the black arrows symbolizes the relative proportion of freshwater discharge from the Mississippi (~66%) and Atchafalaya River (~34%) mouths, with the sum comprising 92% of the total freshwater discharge to the northern GoM.Based on Fofonova et al. (2021).

Fig. 12 .
Fig. 12.(a) Annual and 10-yr anomalies in MSL at Grand Isle, Louisiana.As before, data is missing from 1976 and 1977: the dashed line represents values for those years based on the average of anomalies from Galveston and Pensacola.(b) Annual and 10-yr anomalies in precipitation for the northern Gulf of Mexico drainage basin, including that of the Mississippi River.(c) Annual and 10-yr anomalies in Mississippi River discharge at Vicksburg, Mississippi, above the bifurcation that forms the Atchafalaya River (data from Piecuch et al., 2018).

Fig. 13 .
Fig. 13.(a) Mean annual sea-level pressure anomalies (hPa, colour scale) and 10-m winds (m s − 1 , vectors) for years in the upper quartile of annual anomalies in MSL at Grand Isle (Fig. 9a).(b) Mean annual precipitation anomaly for years in the upper quartile of annual anomalies in MSL at Grand Isle.(c) Mean annual sea-level pressure anomalies and 10-m winds (m s − 1 , vectors) for years in the lower quartile of annual anomalies in MSL at Grand Isle.(d) Mean annual precipitation anomaly for years in the lower quartile of of annual anomalies in MSL at Grand Isle.The Mississippi River Delta study area is shown by the green boxes.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .
Fig. 14.Summary conceptual model for links between the Atlantic Multidecadal Oscillation, 10-yr anomalies in Grand Isle MSL, precipitation in the Gulf of Mexico drainage area, and Mississippi River discharge, and periods of accelerated and reduced land loss.Blue and brown arrows represent years when land-area measurements were made by Couvillion et al. (2017), and which correspond to the AMO negative (cool) or positive (warm) phase, respectively.Yellow arrows represent transitional years.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)