Can Offshore Meteoric Groundwater Generate Mechanical Instabilities in Passive Continental Margins?

Offshore meteoric groundwater (OMG) has long been hypothesized to be a driver of seafloor geomorphic processes in continental margins worldwide. Testing this hypothesis has been challenging because of our limited understanding of the distribution and rate of OMG flow and seepage, and their efficacy as erosive/destabilizing agents. Here we carry out numerical simulations of groundwater flow and slope stability using conceptual models and evolving stratigraphy—for passive siliciclastic and carbonate margin cases—to assess whether OMG and its evolution during a late Quaternary glacial cycle can generate the pore pressures required to trigger mechanical instabilities on the seafloor. Conceptual model results show that mechanical instabilities using OMG flow are most likely to occur in the outer shelf to upper slope, at or shortly before the Last Glacial Maximum sea‐level lowstand. Models with evolving stratigraphy show that OMG flow is a key driver of pore pressure development and instability in the carbonate margin case. In the siliciclastic margin case, OMG flow plays a secondary role in preconditioning the slope to failure. The higher degree of spatial/stratigraphic heterogeneity of carbonate margins, lower shear strengths of their sediments, and limited generation of overpressures by sediment loading may explain the higher susceptibility of carbonate margins, in comparison to siliciclastic margins, to mechanical instability by OMG flow. OMG likely played a more significant role in carbonate margin geomorphology (e.g., Bahamas, Maldives) than currently thought.

Initially ignored because of the wider acceptance of the turbidity current paradigm (Daly, 1936), the groundwater hypothesis was revitalized in the 1980s and 1990s (Orange et al., 1994;Paull et al., 1990;Robb, 1984Robb, , 1990. OMG flow has been proposed to play a role in a range of seafloor geomorphic processes ( Figure 1; Table 1). In siliciclastic margins, OMG may be associated with submarine slope failures by generating excess pore pressures and lowering effective stress (Boffo et al., 2020;Kelner et al., 2016;Kopf et al., 2016;L'Heureux et al., 2010;Oehler et al., 2017;Paull et al., 2021;Stegmann et al., 2011;Sultan et al., 2020). Pore pressure fluctuations and focused fluid flow can lower the consolidation state of sediments and deform them via processes such as creep, liquefaction, fluidization, piping, hydrofracturing, and the development of shear zones (Bull et al., 2009;Kopf et al., 2016;Moernaut et al., 2017;Nardin et al., 1979). In carbonate margins, overpressure development in sub-seafloor aquifers during a relatively rapid sea-level fall has been proposed as a trigger of large-scale failures, which, in low-angle submarine slopes, are thought to have emplaced megabreccias (Busson et al., 2021;Spence & Tucker, 1997). In volcanic islands (e.g., Hawai'i), Iverson (1995) suggested that topographically controlled meteoric groundwater and sea-level fluctuations can only generate the groundwater seepage forces required to cause flank collapse in special circumstances (e.g., the occurrence of a buried clay layer under the growing weight of a volcanic edifice).
Where it actively discharges at the seafloor, OMG has been related to the development of pockmarks and semi-circular depressions (e.g., Andresen et al., 2021;Goff, 2019;Hillman et al., 2015;Hoffmann et al., 2020). Groundwater discharge is thought to exert a seepage force on sediments that offsets the frictional forces within the sediments and causes them to fail or to be eroded grain by grain (Pratson et al., 2009;Robb, 1990). The recurrence of these processes results in a pit. If seepage failure undercuts the pit headwall, the pit can grow upslope by retrogressive slope failure, resulting in a submarine canyon (Green & Uken, 2008;Puig et al., 2017;Robb, 1990). OMG flow and seepage in carbonate margins have been associated with dissolution structures such as sinkholes and conduits, particularly during sea-level lowstands, and other karst landforms such as solution pits, rillenkarren, erosional benches, and tafone recesses (e.g., Back et al., 1984;Land et al., 1995;Lofi et al., 2012;Rousakis et al., 2013). The loss of support at the base of the continental slope due to dissolution of chalk and claystone has been proposed to lead to box canyon development (Orange et al., 1994;Robb, 1984).

Knowledge Gaps
There are two key knowledge gaps that challenge the proposed links between OMG flow and seafloor processes in the above studies.
First, our understanding of the distribution and rate of OMG flow and seepage, as well as their pressure/chemical characteristics, is limited. Information on OMG is predominantly derived from incidental discoveries in boreholes (Micallef, Person, et al., 2021). There is a paucity of direct field measurements of these properties (e.g., Attias et al., 2020;Gustafson et al., 2019;Haroon et al., 2021;Micallef et al., 2020;Sultan et al., 2020), whereas active OMG seepage is difficult to detect and measure (Taniguchi et al., 2019). A cost-effective solution to address these problems is the use of palaeo-hydrogeological models (e.g., Cohen et al., 2010;Micallef et al., 2020;Thomas et al., 2019). Such models provide insights into the mechanisms of OMG emplacement, its distribution and its characteristics. The main issue with published palaeo-hydrogeological models is that they are based on a static grid that does not incorporate the interaction between evolving stratigraphy, hydrology, and sea level. Therefore, we have a poor understanding of how the evolution of a continental margin during multiple sea-level cycles has influenced OMG dynamics.
Second, the efficacy of OMG flow as an erosive/destabilizing agent still needs to be validated. There is a lack of mechanistic understanding of these geomorphic processes, their rates and their spatio-temporal scales. It is still unclear if and what role OMG plays in the formation of submarine landforms, and under which conditions. OMG-driven processes are extremely difficult to measure and observe, and their reproduction with numerical models or laboratory simulations has been rare. All of this inhibits our ability to reconstruct and predict seafloor evolution using groundwater-related processes. during a late Quaternary glacial cycle can generate the pore pressures required to trigger mechanical instabilities on the seafloor. We define a mechanical instability as the mobilization of sediment induced by the effective stress of the material exceeding the yield strength, which can manifest as slope failure and seafloor erosion. Our assessment is based on numerical simulations of groundwater flow and stability carried out with conceptual models and evolving stratigraphy for cases of passive, non-glaciated siliciclastic and carbonate margins. We focus on the latter types of margins because they are the most common types of continental margins globally, where the majority of OMG is hosted (Micallef, Person, et al., 2021), and where the majority of the landforms in Figure 1 is located. Convergent and glaciated margins, in comparison, are influenced by additional geological and hydrological factors that make deciphering the geomorphic role of OMG more difficult.

Significance
Landslides, canyons, and pockmarks are the main drivers of geomorphic change in continental margins. Slope failures are a very common process on siliciclastic margins (Moscardelli & Wood, 2015), mixed siliciclastic-carbonate margins (Puga-Bernabéu et al., 2022), and the flanks of carbonate platforms (Reijmer et al., 2015), including in the geological record (Le Goff et al., 2020;Lehrmann et al., 2020;Spence & Tucker, 1997). They are the key process by which sediments, carbon and nutrients are transported from the shelf to abyssal plains (Masson et al., 2006). Submarine canyons, on the other hand, incise one-fifth of continental margins worldwide (P. T. Harris et al., 2014) and provide a conduit for the transfer of sediments into deep water settings (P. T. Harris & Whiteway, 2011). Sea level has been lower than today for 80% of the Quaternary period (Bintanja et al., 2005). The potential of OMG to influence geomorphic processes across continental margins is likely to have been higher during the majority of the last 2.6 Ma than it is today. Current and past continental margin morphology and processes can  Table 1). Also included are known occurrences of OMG on shelves and slopes (Micallef, Person, et al., 2021), and ice cover during the Last Glacial Maximum (Ehlers & Gibbard, 2003).  Christodoulou et al. (2003) 10.1029/2022JF006954 5 of 25 thus only be understood in terms of a framework where sea levels are lower and OMG systems are more extensive than at present. To develop realistic concepts of continental margin development, it is crucial that these hydrological factors, and the associated geomorphic processes, are taken into consideration.

Conceptual Models
Conceptual numerical models were developed to carry out controlled numerical experiments and assess how evolving OMG characteristics associated with fluctuating sea levels during the last glacial cycle can impact the geomechanical stability of a passive continental margin. In the conceptual models, an idealized continental margin setting was constructed that is representative of typical siliciclastic margins globally. To reflect the properties and features of realistic margins in our numerical scenarios, we controlled the geometry (i.e., characteristic lengths of the shelf and slope, and the depth of the shelf break), and spatial heterogeneity related to the compaction and variation in sediment type from the coast to sea. Uniform and continuous sediment burial were considered over time.
The following model was solved for sub-surface flow fields and the stresses they induce in the sediment matrix: where, ( ) is the sediment porosity, ( ) is the pore-water pressure, , ( , , ) are the water and sediment phase compressibilities, , ( , , ) are water and sediment phase densities, is the acceleration due to gravity, ( ) is the effective Cauchy stress tensor, ( , , ) ∶= − ⊢(∇ + ) is the Darcy velocity of the water phase with permeability ( ) and dynamic viscosity of water , and finally, ( ) ∶= [0 ] , such that: where 0 is the effective burial rate measured at the seafloor. : = denotes the definition of a new variable.
Equation 1 corresponds to the total mass balance of water and sediment phases, Equation 2 to the mass balance of the sediment phase, and Equation 3 to the quasi-steady state momentum balance of the sediment matrix. The Cauchy stress is resolved using an elastic constitutive law in the limit of infinitesimal strains. The coupled system of governing Equations 1-3 is solved for the variables P = [P, ϕ, σ] T . The potential for mechanical failure is measured through the estimation of failure-states (F): Equation 6 denotes the Drucker-Prager yield surface (Φ ) with mean stress ∶= 1 2 ( ) and shear stress ∶= 1 2 ∶ , coefficient of frictional resistance α and coefficient of cohesive strength c, while Φcr denotes the critical yield stress, which is chosen to be a rather conservative value of 10 kPa in this study, based on the average lower bounds for marine sand-mud mixtures (Chen et al., 2021). Failure states less than 0 are considered stable and those above 1 are considered critical. The critical failure states only indicate an expectation for some form of mechanical instability, and the prediction of the actual mode of failure requires further modeling of the underlying geomechanics, which is not considered here.
The computational domain consists of a continental margin with an idealized topography ( Figure 2a) characterized by the length of the shelf L m ∈ {40,90,120}, the length of its slope L s ∈ {14,28,57} km, depth of the shelf break H sb ∈ {70,120,220} m, and the depth of the slope H s = 2 km. Figure 2b) is imposed on the seafloor boundary as Dirichlet pressure constraints. On the exposed areas of the shelf (i.e., areas lying above sea level), a constant average rain rate R (cm/a) is assumed. The corresponding meteoric recharge rate is specified as some percentage r ∈ {1,25,50}% of R.

Sea-level change (shown schematically in
The sediment type is assumed to show a smooth gradient from sandy material coastwards to clayey material seawards. The heterogeneity in corresponding hydraulic and mechanical properties M ∈ {ϕ, κ, E, ν, α, ρ s , C s } is described using the following mapping (shown schematically in Figure 2d): with, In the set , is the evolving sediment porosity, ∶= 0 is the sediment permeability tensor, where − relationship is parameterized using an exponential model (e.g., Hommel et al., 2018;Rutqvist et al., 2002), are the Young's modulus and Poisson's ratio respectively, is the phase density, and is the compressibility of the homogenized granular sediment. Moreover, as the margin forms over time through continuous (or episodic) deposition and burial processes, we consider an anisotropic permeability distribution, rotated along the topography as shown in Figure 2c, to reflect the layered stratigraphy of a typical margin in the ideal limit of continuous burial of identically graded material at a constant rate: where, 0 is the scalar permeability of the sediment, ≤ 1 reflects the degree of flow anisotropy, and ( ) is the angle of rotation (relative to the x-axis) of the margin topography.
In total, 729 scenarios were considered, which include all combinations of the parameters listed in Table 2. The corresponding material properties and model parameters chosen for these simulations are listed in Tables 3 and 4.

Models With Evolving Stratigraphy
In contrast to the conceptual models, the models with evolving stratigraphy consider more realistic representations of two types of passive continental margins-siliciclastic and carbonate-together with associated evolution of groundwater flow and slope stability. The siliciclastic and carbonate margin models were developed using parameters from offshore New Jersey (NE USA) and the western Great Bahama Bank, respectively. The latter is an example of a tropical isolated carbonate margin. The goal of this part of the study was to develop models of stratigraphy and groundwater that are similar to, rather than exact replications of, the field observations from these two study sites.

Data
The parameters for the siliciclastic margin were derived from the following data: (a) bathymetry: Coastal relief model, National Centres for Environmental Information (https://www.ncei.noaa.gov/); (b) seismic reflection profiles: 2-D profiles acquired during expedition OC270 (Mountain, 2008); and (c) boreholes from IODP is the length of the shelf, the length of the slope, and sb the depth of the shelf-break measured with respect to the maximum sea level. The datum for depth z (i.e., z = 0) is chosen to coincide with the depth of the shelf-break. In panel (b), Δ sb = 120 m is the total change in sea level over one glacial cycle, starting at time 0 = 120 ka before present where the sea level was the highest at 0 , decreasing to a level 0 − Δ SL at LGM = 20 ka before present (Last Glacial Maximum, and finally increasing back to 0 in modern time. The sea-level changes are modeled as a linear ramp function over time. In panel (c), and denote the angle of rotation and permeability tensor for the shelf respectively, and and denote the same for the slope. The relationship between (.) − (.) is described through Equation 8. In panel (d), the heterogeneity in the sediment properties due to the sand-to-clay distribution from coast toward offshore is modeled as a quadratic mapping (see Equation 7).

Forward Stratigraphic Modeling
A 3D process-based forward sedimentary-stratigraphic numerical model based on SIMSAFADIM (Bitzer & Salas, 2001, 2002 and SIMSAFADIM-CLASTIC (Clavera-Gispert et al., 2012Gratacós et al., 2009) was used to construct stratigraphic architecture and facies distribution in continental margins. The main sedimentary processes that are modeled by the code are fluid flow, sediment transport, and deposition. The simulation for the siliciclastic margin includes the transport and sedimentation of inflowing terrigenous clastic sediments, whereas for the carbonate margin it also simulates processes of autochthonous marine carbonate production and accumulation by modeling the carbonate producing organisms' evolution and interaction with species associations. The modeling framework, system equations and parameters are presented in detail in Bitzer and Salas (2001), Bitzer and Salas (2002), Gratacós et al. (2009), Gratacós et al. (2021, and Clavera-Gispert et al. (2017). Subsidence and sea-level variation is also taken into consideration (based on Hansen et al., 2013). The code was implemented in FORTRAN 95 and a finite element method was used to discretize the modeled basin and solve the equations of the processes considered.
The required parameters to set up the margin profiles for the stratigraphic modeling were extracted from the integration of stratigraphic interpretation of the available seismic reflection profiles with age and petrophysical Type of sediment at coast --Sand, silt, and clay Hydraulic and mechanical properties of for each sediment type are listed in Table 3 Note. Each combination of the listed parameters constitutes one scenario.  (2001) Note. The subscript "0" denotes the value at the surface (i.e., = 0 ).
information from the exploration well logs. As SIMSAFADIM-CLASTIC works in 3D, an initial modeled area was defined according to a normal-to-shoreline long narrow area (250 × 4 km for the siliciclastic margin, and 33 × 2 km for the carbonate margin) with initial basin topography as a function of the interpreted seismic profiles. The siliciclastic margin profile modeled area was discretized into a finite element mesh with a spatial resolution of 1 × 1 km. The final stratigraphic infill resolution was obtained considering the total modeling time (30 Ma) and the time steps defined (195). The modeled sedimentary units were sand, silt and clay. The carbonate margin profile modeled area was discretized using a spatial resolution of 250 × 500 m. The total modeling time was 1.95 Ma, divided into 195 time steps. This shorter time span was chosen to encapsulate a significant progradation pulse that has been documented for the Bahamian platform (Busson et al., 2019). In this case, the modeled sedimentary units were drift, mass transport deposit, platform, slope, and channel. Sedimentation rates for the siliciclastic margin (Pleistocene: 0.011-0.34 and Oligocene-Miocene-Pliocene: 0.0001-0.02 mm/a) were derived from Dugan and Flemings (2000). For the carbonate margin, sedimentation rates of 0.21-0.46 mm/a for the Pleistocene were derived from Eberli et al. (1997).

Groundwater Modeling
We applied a groundwater flow and solute transport model to represent the hydrogeologic evolution of the siliciclastic and carbonate margins using RIFT2D. The following variable-density groundwater flow equations were solved using a modified version of RIFT2D (Wieck et al., 1995): where S s is specific storage, h is freshwater hydraulic head, ∕ is the total derivative, is the one-dimensional loading efficiency, ∕ is the change in the vertical load through time (Pa/yr) resulting from sedimentation, ⃖ ⃗ is the Darcy flux, K is the hydraulic conductivity tensor, z is elevation, is fluid density, is the relative , is the base density at standard temperature, is the relative fluid viscosity ( = ∕ ), is the base fluid viscosity at standard temperature, concentration, and pressure (10°C, 0 mg/l, atmospheric pressure).
The vertical load is related to the sedimentation rate and effective density of the porous medium: where is the sedimentation (or erosion) rate, g is the gravity constant and is the sediment density.
In order to account for variable-density effects on groundwater flow (Post & Kooi, 2003), the following solute transport equation was solved: where is porosity, C is concentration, and D is the hydrodynamic dispersion-diffusion tensor. The components of the hydrodynamic dispersion-diffusion tensor (in two dimensions: D xx , D zx , D xz , D zz ) are given by Bear (1972): where is the solute molecular diffusivity, and are the longitudinal and transverse dispersivities, and and are the seepage velocities in the x-and z-directions, respectively. The seepage velocities are related to the Darcy flux ( = ∕ ; = ∕ ). Polynomial equations of state (not shown), presented in Batzle and Wang (1992), were used to relate the total dissolved solid concentration to pressure, and temperature to fluid density and viscosity. A conductive-convective heat transfer equation (not shown) was solved, which resulted in a subsurface temperature gradient of about 30°C/ km for use in the equation of state calculations. Convective effects on computed temperatures were minimal.
No flux boundary conditions were imposed on the side and base of the model domain. At the top surface of the model domain, a specified head and concentrations were imposed. If sea-level elevation was above the top nodal elevation of a given column, then the specified head was set equal to sea level and accounting for seawater depth: where bc( ) is the specified head for nodes below sea level, ( ) is the water depth below sea level at a given location and time, and (t) is sea level elevation. If the sea level was below the top node of a given column, then the land surface elevation was specified as the upper boundary condition for groundwater flow. If the sea level was below the top node of a given column, then the concentration was set to 0 PSU. If sea level was above the top node of a given column, then the concentration was set to 35 PSU.
Deviatoric pore pressures were calculated as follows: where Δ P is the deviatoric pressure of the jth node in the ith column and top is the head at the top of the ith column.
RIFT2D was applied to the representations of passive siliciclastic and carbonate margins derived from the stratigraphic models. The triangular finite element grids generated by RIFT2D were constructed along a series of nodal columns. The siliciclastic and carbonate model grids used 120 and 131 nodal columns, respectively. The total number of nodes at the end of each simulation was 8,520 and 9,108, respectively. The maximum number of nodal rows for the siliciclastic and carbonate margin models was 70 and 68, respectively. The distance separating nodal columns was 1,755 and 250 m, respectively. The vertical element widths varied with maximum widths of 15.8 and 3.9 m, respectively.
Hydrogeologic parameters used for the siliciclastic and carbonate models are listed in Table 5. The hydraulic conductivity values for the siliciclastic margin are consistent with those of Dugan and Flemings (2000) and Cohen et al. (2010), whereas the values for the carbonate margin are consistent with those of Jones et al. (2002). For solute transport, longitudinal and transverse dispersivities of 10 and 1 m, respectively, were assigned to all units.
The duration of the model runs for the siliciclastic and carbonate margins, and the sea-level curves, were the same as the forward stratigraphic model. The simulation time for the siliciclastic margin comprised 30,000 time steps of 1,000 each, whereas that for the carbonate margin consisted of 390 time steps of 5,000 each.
The primary simulations entailed consideration of sea-level fluctuations, sedimentation and mechanical loading for a margin with evolving stratigraphy. An additional two simulations were run. In the first one, sea-level fluctuations were considered but there was no mechanical loading due to sea-level changes or sedimentation. In the second, a static grid scenario was considered where the present-day grid was used and boundary conditions identical to those in the evolving model were imposed. This required an initial salinity distribution to be assumed. The salinity of all nodes below a seawater elevation of −40 m was set to 35 PSU. This represents the average sea-level elevation for the Pleistocene (Meisler et al., 1984). All nodal columns that had an elevation above −40 m were initially assigned 0 PSU. The static grid runs had no subsidence included (and hence, no overpressure formation due to sediment loading).

Slope Stability Modeling
The general limit equilibrium code Slide2 (Fredlund & Krahn, 1977;Rocscience, 2020) was used to estimate the factor of safety (FS) of the continental slopes of the simulated margins (at present and the LGM), and the failure geometry. The equation that was used to calculate the FS is derived from Dugan and Flemings (2002): where θ is the slope angle, λ* is the overpressure ratio (λ* = Δu/δ'vh), u is the pore water pressure (u = ρ w gz), δ'vh is the hydrostatic vertical effective stress (δ'vh = ρsgz), c' is the effective cohesion, φ' is the effective soil friction angle, γ' is the effective soil unit weight (γ' = ρ s g), ρ s is sediment density, ρ w is fluid density, z is sub-seafloor depth, and g is gravitational acceleration. FS ˃ 1 is indicative of stability whereas FS < 1 denotes instability.
The morphological information (e.g., slope angle, depth) was derived from the relevant output of the forward stratigraphic model. The values for the geotechnical properties used for the sediments in the siliciclastic and carbonate margins are listed in Tables 6 and 7. These tables take into consideration the impact of consolidation on the different geotechnical properties. The siliciclastic and carbonate margins were divided into eight and four zones, respectively, on the basis of changes in porosity with depth. Changes in porosity, unit weight and shear strength were directly derived from IODP and ODP boreholes from offshore of New Jersey and the Great Bahama Bank. The friction angle was estimated using the Mohr-Coulomb failure criterion (Terzaghi, 1943): where τ is the shear strength and ′ is the effective normal stress. The values of cohesion for the siliciclastic margin are taken from Table 3, whereas the carbonate sediments are considered to be cohesionless (Lavoie, 1988). Note.
where K x = horizontal hydraulic conductivity; K z = vertical hydraulic conductivity; S s = specific storage.

Conceptual Models
We first report how the various parameters considered in the conceptual models influence the mechanical stability of a continental margin (Figure 3). For this, we have chosen as a reference case the scenario with mid-value for each parameter in Table 2, that is, = 90 , = 28 km, sb = 120 m, = 25% , 0 = 0.025 cm/a, and sand-to-clay material distribution. The trends described below are with respect to this reference case.
1. Impact of length of shelf : Under the assumption of a smoothly varying sediment type, longer margins offer lower gradients of material properties . At the same time, they also offer a larger recharge area, leading to competing feedback on the mechanical stability. For sufficiently high gradients of material properties (e.g., sand or silts at the coast), longer shelves lead to larger regions of mechanical instability (Figure 3b).   2. Impact of length of slope : Steeper slopes (i.e., smaller L s ) provide higher gradients of material properties in the vicinity of the shelf-break as well as larger topographic gradients for fluid flow, and therefore, lead to larger regions of mechanical instability. 3. Impact of depth of shelf break sb : The depth of the shelf-break controls the rate and extent to which the shelf is exposed to meteoric recharge as the sea level drops. For margins where the shelf-break lies close to the lowest sea level, the critical states are expected to appear at the sea-level minima, depending on the gradients of change in granular material properties and the recharge rates. For shallower shelf-breaks, a larger shelf area is exposed as the shoreline recedes faster, leading to the appearance of critical states already before the sea level minima (Figure 3d with sb = 70 m where the shelf-break lies 50 m above the lowest sea level, and the critical states appear nearly 40 ka before LGM). On the contrary, for deeper shelf-breaks, the shoreline recedes at a slower rate and less shelf area is exposed, leading to a lower likelihood of mechanical instability ( Figure 3d with sb = 220 m where the shelf-break lies 100 m below the lowest sea level, and the critical states do not appear). 4. Impact of recharge rate : Higher meteoric recharge directly affects the pressure gradient in the subsurface, and therefore, leads to higher mechanical instability, given sufficiently high gradient of change in granular material properties and exposed shelf area available for recharge ( Figure 3e). 5. Impact of burial rate 0 : Sediment burial leads to compaction (see porosity evolution in Figure S3f in Supporting Information S1), which tends to increase the mechanical stresses. Higher burial rates lead to larger regions of mechanical instability (Figure 3f). Within the conceptual model, we only account for a continuous sediment burial and choose an average value of 0.1 cm/a based on the global burial rate estimates of Wallmann et al. (2012). If sediment loading is large enough, or episodic, we can expect the impacts of burial-related sediment loading to become more dominant than those of sea-level changes. 6. The impact of sediment variability: The gradient of material properties has a direct correlation with mechanical instability. In our simulations, clayey sediments did not show any critical states, even for the combination of longest margin with steepest slope and highest recharge rate considered. Sediments with sand-to-clay variation showed the most margin configurations that could become mechanically unstable during sea-level changes.
The location and timing of mechanical instability are controlled by topography (i.e., a combination of L m , L s , and H sb ) and sediment variability. Within the constraints of the explored parameter space, our model suggests that the highest chance of failure is in the vicinity of the shelf break, at or shortly before the LGM (Figure 3). Broadly, this is because the highest sediment variability can be expected in the vicinity of the shelf-break, with the highest pressure gradients expected around the LGM. We expect mechanical instability before the LGM, rather than after, because our model suggests that the failure states are induced by a decrease in sea level. Figures S1-S4 in Supporting Information S1 display snapshots of pressure distribution, pressure gradients and porosity distribution for selected scenario margin configurations. Figure 4a presents the configuration and hydraulic conductivity distribution for the siliciclastic margin model at present-day. The hydraulic conductivity decreases in an offshore direction, with the highest hydraulic conductivities recorded in the upper shelf, which was subaerially exposed during most of the Pleistocene. The vertical alternation of coarse sand and silt facies can be observed in the mid-shelf, and is the result of sea-level fluctuations. Figure 4b presents the configuration and hydraulic conductivity distribution for the carbonate margin at present-day. The highest hydraulic conductivities are found at the top of the carbonate platform. Hydraulic conductivity decreases with depth on the platform, and laterally along the platform margin and slope facies. The alternation of coarse and fine facies occurs between the shelf break and the slope, with the least permeable facies dominating the deep-water deposits.

Groundwater Flow Models
Computed deviatoric fluid pressures for the siliciclastic margin at the LGM and present-day are displayed in Figures 5a and 5b. There is very little difference between the simulated heads at the LGM and the modern conditions. The maximum heads occur near the base of the model domain where the low permeability offshore deposits are thickest. The maximum deviatoric pressure is about 10 MPa. Figures 5c and 5d present computed deviatoric pressures without sediment or sea-level loading for the LGM and present-day. The deviatoric pressures are considerably lower when sediment and sea-level loading is neglected. The highest computed deviatoric pressure without sediment loading is on the order of 0.8 MPa. This is controlled by topographically driven flow across the continental shelf.
Computed deviatoric pressures for the carbonate margin are significantly lower than those of the siliciclastic margin (Figures 5e and 5f). The top of the carbonate platform is located at a depth of 38 m. During the LGM, a recharge area formed with descending groundwater flow within the carbonate platform that produced negative deviatoric pressures (Figure 5e). Positive deviatoric pressures occurred below the shelf break near the intersection of the break deposits and sea level. At present-day, the carbonate platform is flooded, and changes in head due to variable-density flow control the distribution of deviatoric pressures (Figure 5f). When we removed the sediment and sea-level loading terms from the model, the resulting computed deviatoric pore pressures were nearly identical (results not shown).
Computed salinity patterns for the siliciclastic and carbonate margin simulations are presented in Figure 6. In the evolving siliciclastic margin models, the greatest endowments of OMG occurred during the LGM (Figure 6a), when a wide expanse of the continental shelf was exposed to meteoric recharge. Freshened groundwater (<33 PSU) extended out to the continental shelf break at a sub-seafloor depth of 50-500 m during the LGM. As sea level rose to present day conditions, the groundwater near the sediment-water interface was salinized within the permeable sand and silt facies due to haline convection and diffusion (Figure 6b). The penetration of saline water was <100 m. Due to the long response times for solute transport, the freshened groundwater did not change much since LGM times.
Most models of continental margin hydrogeology do not represent variations in sedimentation and subsidence during the Pleistocene. Figure 6c shows the result obtained from an additional scenario using a static grid and no subsidence, where initial conditions were assigned. In contrast, the initial salinity condition for a given element in the evolving grid model (Figure 6b) is generated at the time when the element is created at the sediment-water interface. The local salinity conditions control what salinity level is assigned to each element. We observe differences in the present-day salinity patterns between the evolving grid and static grid model results (Figures 6b  and 6c). In particular, for the evolving grid, there is more saline water at depth that has not been flushed during the Pleistocene. This is due to the relatively slow response time for advection and vertical diffusion/dispersion at great depths.
For the carbonate margin simulations, the groundwater in the upper permeable facies of the platform was entirely fresh during the LGM (Figure 6d). As a result of sea level rise that led to present-day conditions, vertical diffusion and haline convection salinized the upper platform, especially within the permeable upper 30 m (Figure 6e). Freshened groundwater is preserved in the permeable platform facies at depth. Sub-seafloor at a water depth of >120 m is dominated by high salinity conditions. The static grid model result does not appear to be significantly different from that of the evolving grid model (Figure 6f). This may be because the model domain is relatively thin when compared to the siliciclastic margin.

Slope Stability Models
The estimated FS for the siliciclastic margin is >1 for present-day and LGM conditions for scenarios with and without sediment and sea-level loading (Figures 7a-7d), indicating continental margin stability. The computed FS is lower during LGM than under present-day conditions, and higher for the simulations ignoring sediment and sea-level loading. For the carbonate margin, the FS is >1 for present day conditions, but decreases to 0.9 during the LGM, suggesting instability (Figures 7e and 7f). The section of the margin that is estimated to have been susceptible to failure is 1 km wide and located in the upper slope.

Comparison With Field Data
Although our models with evolving stratigraphy were not meant to replicate the continental margins of New Jersey and the Great Bahama Bank, we do observe general similarities in terms of stratigraphy, groundwater characteristics, and pore pressures.
For the siliciclastic margin, the hydraulic conductivity distribution in Figure 4a is comparable to the hydrostratigraphic model developed by Meisler et al. (1984) for New Jersey, with hydraulic conductivity decreasing overall in an offshore direction. The salinity distribution across the New Jersey margin has been represented as a salinity contour map in Meisler et al. (1984), which was generated from borehole salinity profiles in Hathaway et al. (1979), and the salinity profiles in Lofi et al. (2013) (Figure 6b). In these representations, a laterally extensive fresh-brackish water lens is located at a sub-seafloor depth of 50-400 m, with salinity increasing at larger depths. The salinity distribution in Figure 6b is generally comparable to these representations, although we do observe an underprediction of the OMG thickness in the middle shelf and an overprediction in the outer shelf. Computed overpressures are similar to those simulated by Dugan and Flemings (2000). This is encouraging, considering that we used the same compressibility and vertical permeability as reported by these authors.
For the carbonate margin, the hydraulic conductivity distribution in Figure 4b is similar to the permeability distribution modeled for the Great Bahamas Bank by Busson et al. (2021). We do not have information on groundwater distribution on the platform, but on the slope the modeled salinity distribution in Figure 6e matches salinity measurements in boreholes from ODP Expedition 166.

OMG Flow as a Driver of Mechanical Instability
In our simulations, enhanced topographically driven flow during sea-level lowstands, resulting from high hydraulic heads and increased recharge across an exposed continental shelf, generated higher pore pressures than at LGM and (f) present-day with sediment and sea-level loading. The extent of margin displayed in these figures is shown in Figure 4. The red curved line denotes the critical slope surface. present in the outer shelf to upper slope sediments. These pore pressures predominantly developed in zones of high sediment variability and reduced the materials' shear strength, giving rise to instability. Both the conceptual and evolving stratigraphic model results (Figures 3 and 7) show that mechanical instability by OMG flow is most likely to occur in the outer shelf to upper slope region, having an extent of 1-6 km, at or shortly before the LGM sea-level lowstand. The shallower the shelf break depth, the earlier the timing of the instability is. For models with evolving stratigraphy, which are more realistic representations of continental margins than the conceptual models, it is only in the carbonate margin that OMG flow is the main driver of pore pressure increase and associated instability. In the siliciclastic margin, the main process driving an increase in pore pressure and preconditioning the slope to failure is sediment loading (in agreement with Dugan and Flemings (2000) for the New Jersey setting), with OMG flow playing a comparatively minor role. Potential reasons for this difference are explored in Section 4.3.

Factors Controlling Mechanical Instability by OMG Flow
There are three main factors that can explain the higher susceptibility of passive carbonate margins to instability by OMG flow: 1. Sediment variability: The conceptual model results ( Figure 3) indicate that mechanical instability by OMG flow in passive margins is promoted by: (a) high recharge rate, which tends to increase with the length of shelf and precipitation rates and (b) high sediment variability, which tends to be inversely related to shelf/ slope length. The relevance of sediment variability is also evident in the results from models with evolving stratigraphy (Figure 7). Carbonate margins, especially those dominated by autochthonous sediments, have a stronger degree of spatial and stratigraphic heterogeneity in comparison to siliciclastic margins (Dalrymple & James, 2010;Tucker & Wright, 1990). Carbonate margins comprise sediments ranging from boulders to mud. Siliciclastic margins, in comparison, comprise a narrower range of grain sizes (sand to mud) and an increase in low-permeability sediments from nearshore to offshore. 2. Sediment strength: The shear strength of carbonate sediment in our models tends to be lower than that of siliciclastic sediments (Tables 6 and 7). This is in line with the general trend observed in global scientific drilling sample measurements (Bartetzko & Kopf, 2007), and would explain the higher susceptibility of carbonate margins to failure. 3. Overpressure development by sediment loading: Our results suggest that there is a competition between sediment loading and OMG in increasing pore pressures in continental slopes. Overpressures due to sediment loading tend to be higher on the slopes of siliciclastic margins because: (a) sedimentation rates are generally higher than those in carbonate slopes (Einsele, 2000) and (b) siliciclastic sediments tend to have high compressibilities, whereas compressibility and overpressure development during carbonate sediment burial is generally limited by cementation and biological binding (Armitage et al., 2018;Dugan & Flemings, 2000;Moshier, 1989).
Although not considered in our models, there may be two additional reasons why carbonate margins are more susceptible to mechanical instability by OMG flow. First, water-rock interactions in carbonate margins (e.g., dissolution, micritisation, dolomitization) can give rise to multi-scale, heterogenous permeability structures that may enhance the extension of groundwater offshore, even under present conditions (e.g., Evans & Lizarralde, 2003;Legrand & Stringfield, 1971;Sanford & Konikow, 1989;Swarzenski et al., 2001). Second, cementation and biological binding in carbonate margins can increase the potential for brittle failure and gravitational collapse (Armitage et al., 2018;Moshier, 1989). In siliciclastic margins, early lithification is uncommon, and the potential for brittle failure is low.

Global Significance of OMG Flow as a Geomorphic Agent
The seafloor landforms attributed to OMG in Table 1 predominantly occur in siliciclastic margins in shallow settings close to the coast. This is a different setting from that modeled in our study, and active groundwater recharge at present via artesian confined aquifers that extend below submarine regions may develop the pore pressures required to generate mechanical instabilities (Sultan et al., 2020). On siliciclastic continental slopes, the role of OMG is likely to be less important in forming pockmarks and canyons than initially suggested (Green & Uken, 2008;Hillman et al., 2015;Hübscher & Borowski, 2006;Puig et al., 2017;Rise et al., 1999), unless elements that promote the extension of OMG toward the outer shelf and its influence on the mechanical instability are/were present. These include well-connected high permeability zones, preferential flow pathways (e.g., buried palaeochannels) or sub-seafloor aquifers exposed by erosion (Groen et al., 2000;Houben et al., 2018;Lin et al., 2010;Michael et al., 2016;Paldor et al., 2020). Offshore New Jersey, regularly spaced blind canyons have been documented in gently seaward dipping chalky Eocene sedimentary rocks and silty Miocene claystones (Orange et al., 1994;Twichell & Roberts, 1982). The formation of these canyons has been attributed to groundwater seepage during sea-level lowstands (Orange et al., 1994;Robb, 1984). We do consider groundwater (not necessarily OMG) to have played a key role in the formation of these blind canyons. However, in view of their occurrence at the base of the continental slope and the observation of karst-like features in calcareous lithologies, it is likely that the process responsible for their formation entailed a reduction of rock strength via fluid pressure and dissolution, as proposed by Micallef, Paull et al. (2021) for the development of box canyons in carbonate escarpments.
The landforms associated with OMG in carbonate settings are primarily linked to dissolution (Table 1), which was not addressed in our study. However, it is in these settings that OMG is likely a major player in creating mechanical instabilities. The heads of blind submarine canyons in the Little Bahama Bank, for example, are located at water depths of ∼450 m, and these are thought to have been initiated by slope failures associated with sediment overloading (Mulder, Ducassou, Gillet, et al., 2012;Recouvreur et al., 2021;Tournadour et al., 2017). Several examples of the latter are observed along the middle slope of the Little Bahama Bank, with pockmarks located above the slide scars (Principaud et al., 2015). Upper-to mid-slope failures have been reported in the NW and SW corners of the Great Bahama Bank (Jo et al., 2015;Principaud et al., 2017) and are thought to have been triggered during eustatic falls after periods of very high rates of slope sedimentation. Here, pockmarks have documented upslope of the scars (Mulder, Ducassou, Eberli, et al., 2012). Slope failures along the slopes of the Maldives carbonate platform, on the other hand, have been attributed to a reduction in sediment shear strength moderated by fluid ascent and the destabilization of gas hydrates during sea-level lowstands (Lüdmann et al., 2022). Based on our results, we suggest that OMG could also have been a key factor in the development of the landforms documented in uncemented slopes offshore the Bahamas and the Maldives during sea-level lowstands.

Temporal Distribution of Slope Failures
A number of studies have proposed that slope instability, particularly in the Atlantic Ocean, is more frequent during periods of glaciation and/or during glacial to interglacial transitions (Hilbrecht, 1989;Lee, 2009;Maslin et al., 2004;Owen et al., 2007;Paull et al., 1996). Maslin et al. (2004) and Paull et al. (1996) suggest that this is mainly associated with lowering sea levels because of reduced hydrostatic pressure and the associated destabilization of gas hydrates. Lee (2009), on the other hand, links the higher frequency of slope instability to the development of thick sedimentary deposits on the upper continental slope during sea-level lowstands, which fail due to seismicity triggered by isostatic readjustment the of previously glaciated regions.
In the Mediterranean Sea, which predominantly hosts active margins, most of the largest magnitude failures were mobilized during sea-level lowstands (Urgeles & Camerlenghi, 2013). These failures were associated with pore pressure generation driven by high sedimentation rates and gas exsolution. Similar patterns have been observed in turbidity current activity offshore Portugal, Mauritania, and the Arabian Sea (Henrich et al., 2010;Lebreiro et al., 2009;Prins et al., 2000). The higher likelihood (3-5 times) of slope failure during falling or lowering sea level, in comparison to rising or high sea level, is supported by modeling results (Hutton & Syvitski, 2004). Urlaub et al. (2013), on the other hand, used a global compilation of 41 large submarine landslide ages in the last 30 ka to document a weak global correlation of landslide frequency with sea-level changes or increases in local sedimentation rate, which they attributed to uncertainties in the landslide age. Pope et al. (2015) conclude that there are currently too few, sufficiently well-dated large landslides to determine whether they are temporally random.
Based on our results we surmise that more frequent slope instability during sea-level lowstands, if proven to be statistical significant, could also be indirectly attributed to OMG flow in certain settings.

Limitations and Outstanding Knowledge Gaps
The scope of our study was limited to exploring whether OMG flow can generate mechanical instabilities in siliciclastic and carbonate margins. There are a number of limitations associated with our modeling approach: 1. One main difficulty is determining appropriate boundary conditions considering the long time scales being simulated and the general lack of field data. For this reason, the onshore component of the groundwater systems was not considered, and we were unable to take into account the variability of climatic parameters (e.g., precipitation, temperature, land use, and ice cover). 2. The linear elements in our groundwater model do account for sediment anisotropy. However, due to the skewness of the elements, errors may be introduced (e.g., Balan et al., 2022;Huang et al., 2017;Kopteva, 2020aKopteva, , 2020b. Additional numerical simulations we carried out (Figures S5-S7 in Supporting Information S1) suggest that there is no evidence of numerical errors based on grid refinement. 3. The generated stratigraphic models are simplified and do not include the fine-scale heterogeneity observed in field data. The models use a single permeability value for each facies, and, being 2D, do not take into account variability in shore-parallel sediment and flow regimes. They also ignore cementation, although this was shown to have a second-order influence on the pore fluid pressure regime and instability in carbonate settings (Busson et al., 2021). 4. The critical yield strength changes with depth, but we have ignored this effect in the slope stability models.
Typically, the yield strength will become higher with depth, and the failure envelope will therefore expand as well (i.e., failure will occur at higher stresses). Since, in our calculations, the failure states are confined to the shallower sediment layers, we expect that ignoring the depth dependence of critical yield stress will not change the results significantly. Moreover, the depth-dependent parameterization would have introduced additional fitting parameters and therefore more uncertainty. 5. We have not considered the impacts of onshore and offshore terrain changes caused by sedimentation and erosion, and rather focused only on sea-level changes. Feedback between terrain change and groundwater dynamics is complex and will be considered in future studies. 6. The Bahamas are an example of a tropical isolated carbonate margin. The outcomes of our study are therefore not directly transferable to other types of carbonate margins (e.g., cool-water carbonate margins such as the Great Australian Bight), or mixed carbonate-siliciclastic margins (e.g., NE Australia).
A number of knowledge gaps remain to better constrain the efficacy of OMG as a seafloor geomorphic agent on a global scale. These include (a) assessing the role of dissolution in carbonates and salt leaching in siliciclasts in changing the permeability architectures of margins and triggering slope failure; (b) determining the role of OMG in driving slope failure in convergent margins (e.g., R. N. Harris et al., 2013), glaciated margins (e.g., Paull et al., 2021), and volcanic islands (e.g., Iverson, 1995); glaciated margins are of particular interest because estimates of OMG recharge rates can be 2-10 times higher than modern values (Person et al., 2007); (c) exploring how OMG flow develops submarine landforms in 3D, and how these landforms may in turn influence OMG flow dynamics; and (d) constraining the style of slope failure associated with OMG flow. A strategy that integrates laboratory experiments with a numerical approach that builds on that used in our study (e.g., by integrating fault and geochemical modeling) and that is extended to 3D, would be suitable to address these knowledge gaps.

Conclusions
We have carried out numerical simulations of groundwater flow and slope stability using conceptual models and evolving stratigraphy to assess whether OMG and its evolution during a late Quaternary glacial cycle can generate the pore pressures required to trigger mechanical instabilities in the passive, non-glaciated siliciclastic and carbonate margins. Our main conclusions are the following: 1. Mechanical instability by OMG flow is most likely to occur in the outer shelf to the upper slope region of continental margins, having an extent of 1-6 km, at or shortly after the LGM sea-level lowstand. 2. OMG flow is a key driver of pore pressure development and instability in carbonate margins, whereas it only plays a minor role in preconditioning the slope to failure in siliciclastic margins. 3. The higher susceptibility of carbonate margins to mechanical instability by OMG flow in comparison to siliciclastic margins is attributed to their higher degree of spatial/stratigraphic heterogeneity, lower shear strengths of their sediments, and limited generation of overpressures by sediment loading.
4. OMG is likely to have played a more significant role in carbonate margin geomorphology than currently thought, and it may explain the formation of slope failure scars, blind canyons, and pockmarks during sea-level lowstands in places such as the Bahamas and the Maldives. The role of OMG in forming pockmarks and canyons in siliciclastic slopes is likely to be less significant than suggested in the literature (Table 1), unless these were connected by high permeability pathways to the coast/inner shelf.

Data Availability Statement
Models: The conceptual model code is available at https://doi.org/10.5281/zenodo.7094202. SIMSAFADIM is described in Bitzer and Salas (2001) and Bitzer and Salas (2002), and SIMSAFADIM-CLASTIC in Gratacós et al. (2009), Clavera-Gispert et al. (2012, and Clavera-Gispert et al. (2017), and available from Òscar Gratacós (ogratacos@ub.edu). RIFT2D is described in Wieck et al. (1995) and in the text, and is available from Mark Person (mark.person@nmt.edu The sources of the data used in the various models are provided in the text. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant 677898 (MARCAN)) and the U.S. National Science Foundation (NSF FRES 1925974). We thank Roger Clavera-Gispert for his assistance with SIMSAFADIM and Or Bialik for the useful discussions. We are grateful to Amy East, Angel Puga Bernabéu, and Qiliang Sun for their insightful comments. Open Access funding enabled and organized by Projekt DEAL.