Sensitivity of Subsurface Permeability in Coastal Deltas to Their Morphodynamic and Geomorphic Characteristics

The amount of fresh water moving through coastal deltas worldwide is controlled by the complex subsurface structures within a delta. Morphodynamic influences produced by the feeding river, waves, and tides, in addition to sea level transgressions and regressions, have resulted in deltaic aquifers that are highly heterogeneous. We use 171 unique two‐dimensional morphodynamic models to explore the range of subsurface permeability, hydraulic gradient, and groundwater flux within three end‐member delta types (fluvial, wave, and tidal). We quantify the connectiveness of the subsurface permeability and estimate the horizontal heterogeneity and anisotropy of the permeability. A distance‐based generalized sensitivity analysis is used to investigate the impact morphodynamic influences (fluvial, wave, and tidal), basin conditions (sediment concentration and bathymetric gradient), and geomorphic characteristics (number of channels, shape of the delta plain, and shoreline rugosity) have on the subsurface permeability, hydraulic gradient, and connectivity. We find that the median permeability in deltaic landforms is 4.0 × 10−12 m2 (relating to a hydraulic conductivity of 2.1 × 10−5 m/s), the average hydraulic gradient is 3.9 × 10−4, and the mean specific discharge is 1.3 × 10−8 m/s. High permeability bodies are highly connected and are associated with channelization. Subsurface permeability, hydraulic gradient, and the connectiveness of high permeability areas are sensitive to morphodynamic influences (fluvial, wave, and tidal) and the geomorphic characteristics (number of channels and shoreline rugosity) within a delta. Since morphodynamic influences and geomorphic characteristics are easily identified by looking at the surface of the delta, we suggest that the deltaic subsurface can be characterized by identifying features on the delta surface.

It is well-documented that the morphology of a delta is dependent on three primary morphodynamic influences controlling the deposition of sediment: (a) the fluvial discharge of the sediment-bearing river, (b) wave energy, and (c) tidal energy.Historically, the geomorphology of a delta has been described qualitatively, comparing characteristics of the delta plain and river network to distinguish end-member morphologies (Figure 1; Galloway, 1975).Fluvial deltas, like the Mississippi, are elongated and typically only have a few distributary channels.Wave deltas, like the Sao Francisco, have a smooth shoreline with one main channel cutting the delta plain.Tidal deltas, like the Mekong, have an abundance of sinuous tidal channels that are wider at the mouth and are separated by narrow inland structures.Other basin conditions, including regional bathymetry (Caldwell et al., 2019), the amount and particle size of sediment entering a delta (Caldwell & Edmonds, 2014;Rossi et al., 2016; J. P. M. Syvitski et al., 2010), and antecedent stratigraphy (Geleynse et al., 2011) explain geomorphic differences in deltas with similar proportions of each morphodynamic influence.Recently, geomorphic characteristics of deltas have been quantified, including the shape of the delta plain, elongation and skewness of the plain, the number of distributary channels, shoreline rugosity, channel sinuosity, and channel shape (Caldwell & Edmonds, 2014; J. P. M. Syvitski & Saito, 2007;L. D. Wright & Coleman, 1973).Research linking surface features in a delta, such as the river network, to the subsurface permeability has only considered fluvial deltaic environments (Hariharan et al., 2021;Kolker et al., 2013;Steel et al., 2022;Xu et al., 2021).
We explore how the morphodynamic influences, basin conditions, and geomorphic characteristics of a delta impact the shallow groundwater flux (or specific discharge) through the delta.We use 171 unique two-dimensional models to characterize the geomorphic characteristics (number of channels, shape of the delta plain, and shoreline rugosity) and hydrogeologic conditions (permeability and hydraulic gradient) for the three end-member delta types: fluvial, wave, and tidal.We investigate the range of permeability and hydraulic gradient in each delta type to estimate the groundwater specific discharge through the deltas.We further look at delta permeability by investigating the horizontal heterogeneity and anisotropy as well as the connectivity of high permeability bodies.We use a distance-based generalized sensitivity analysis (DGSA) to investigate the sensitivity of permeability, hydraulic gradient, and connectivity of the high permeability bodies to changes in the morphodynamic influences, basin conditions, and geomorphic characteristics.

Numerical Morphodynamic Modeling
We simulate the formation of the top Holocene layer in modern coastal deltas (formed within the past 8 ky of marine regression) using Delft3D, an open-source, physically based, numerical morphodynamic modeling software (Deltares, 2013).The morphodynamic models are 2D and use a depth-averaged approach to solving the flow and sediment transport equations (Deltares, 2013).We used depth-averaged flow to reduce model runtimes  and because initial results suggest that the horizontal flow field does not have significant variation in the vertical direction.Using a 2D depth-averaged approach has been used in several other delta modeling studies (Broaddus et al., 2022;Caldwell & Edmonds, 2014).Additionally, using depth-averaged morphodynamic models is applicable to deltaic environments with marine influences including waves and tides.The initial model setup is based on the modeling protocol described in Caldwell and Edmonds (2014).The models mimic a sediment-carrying river debouching into an oceanic basin that increases in depth away from the initial shoreline (Figure S1 in Supporting Information S1).The initial shoreline, located along the southern boundary of the model, consists of an erodible beach 5 m above sea level (m asl) and protrudes 5 km into the basin.The beach and the basin floor are 0.5 m thick and consist of homogenously mixed sediments that have the same distribution of sizes as the incoming sediment.A channel is located halfway along the southern boundary and cuts through the initial shoreline.Bathymetry in the basin is shallowest at −5 m asl in the channel at the southern boundary and deepens to the north with a uniform gradient (this bathymetric gradient is changed in each model).The morphodynamic modeling domain is 80 km × 100 km.The simulation represents approximately 50 years of morphological growth (using a morphological updating factor of 200), with time steps every 12 s (0.2 min).This simulation period is able to generate morphodynamic features (e.g., area) similar to what is observed in real-world large deltas.
Boundary conditions in the model are stable over time and consist of a total discharge boundary located at the river channel, an open ocean (Dirichlet-type) boundary along the northern edge of the model, and zeroflux (Neumann-type) boundary along the eastern and western edges and along the bottom layer (Figure S1 in Supporting Information S1).We created 171 unique models by varying the proportion of wave height (Wa) and tidal amplitude (Ti) at the ocean boundary, in addition to varying fluvial discharge (Qav) entering the model through the discharge boundary (Anderson, 2023).Furthermore, three basin conditions are also varied; these include the concentration of sediment entering the model through the stream boundary (Cs), the initial slope of the model representing regional bathymetric gradient (Dgrd), and the median grain size of sediment (Dmm) entering the model (Anderson, 2023).Variations in Qav, Wa, Ti, Cs, Dgrd, and Dmm are based on data collected from 51 deltas around the world and are representative of conditions within the world's largest coastal deltas (Anderson, 2023).The 51 real deltas from which data were collected are not representative of all deltas around the world due to an intentional sampling biased of large deltas; that is, the data set excludes small deltas that are typically wave dominated and that are often little more than a small extrusion from the shoreline (Nienhuis et al., 2020).The input parameters for all 171 models are included in the Data Set S1.
Each of the 171 morphodynamic delta models is categorized as a fluvial, wave, or tidal by the proportion of the morphodynamic influences (model inputs) present.In fluvial models, more than 50% of the influence is derived from the fluvial discharge (Qav) of the incoming stream (Anderson, 2023).In wave models, more than 40% of the influences are generated through wave height (Wa) at the northern basin boundary.In tidal models, more than 35% of the influences are generated through tidal modulation of the water level at the northern boundary; the degree of modulation is referred to herein as the tidal amplitude (Ti).The percentages that define each category were chosen to produce a similar number of models in each of the three categories.Of the 171 models, 59 are fluvial-influenced (referred to herein as fluvial deltas), 55 are wave-influenced (referred to as wave deltas), and 59 are tidal-influenced (referred to as tidal deltas).The remaining 35 models have a combination of influences that does not allow for categorization in any of the end-member groups.

Shape of the Delta Plain
Delta shape (Sh) is a metric describing the elongation of the delta into the receiving basin.The shape of a delta is found by first determining the extent of the delta plain contained within the shoreline.The shoreline is defined using the Open Angle Method (Shaw et al., 2008).We implement the Open Angle Method by defining any cell where the elevation of sediment exceeds the elevation of water as a land cell and any cell where the water elevation is higher than the sediment elevation as a water cell.The Open Angle Method allows for a definition of the shoreline across each river intersecting the shoreline by defining the angle that spans the view into open water; a small open angle will define a shoreline that closely follows that land-water interface while a large open angle will produce a shoreline that is smoother (Shaw et al., 2008).Several open angles were tested; an angle of 45° was chosen for this study because it can accommodate the widest river mouths produced in the models while still preventing over-smoothing of the shoreline.Shape (Sh) is calculated by dividing the maximum delta width (between the eastern and western boundaries) and dividing it by two times the maximum delta length (between the southern and northern boundaries) (Caldwell & Edmonds, 2014): (1) Deltas with an elongated delta plain will have a smaller Sh, whereas deltas that stretch along the shoreline and do not protrude into the basin will have a larger Sh.

Rugosity of the Shoreline
Shoreline rugosity (Sr) is a sinuosity measurement calculated by dividing the length of the mapped shoreline by the length of a smoothed shoreline (Caldwell & Edmonds, 2014).The smoothed shoreline is calculated using a Gaussian filter on the mapped shoreline.A smoothing window of 30 pixels is used as it reduces the sinuosity of the mapped shoreline while still following the overall curve of the delta plain.

Channel Network
The channel network is mapped using Rivamap (Isikdogan et al., 2017), which identifies individual river segments using a multiscale singularity index to extract curvilinear structures from the elevation map of the delta.
The number of channels (Cn) is calculated by identifying all the river segments that are longer than one-tenth the length of the delta plain.This minimum length ensures that single cells or small depressions in the elevation map are not identified as a river.

Hydraulic Conductivity and Permeability
Hydraulic conductivity (K) of the sediment deposited in the delta is calculated using the Hazen method (Hazen, 1911)  ( where K ij is the hydraulic conductivity at a cell with index [i, j], C is a shape factor (constant), and d 10 is the grain size in which 10% of the sediment is finer.A shape factor of 40 is used in this study, representing the bottom of the range for very fine, poorly sorted sand (Hazen, 1911).This value is selected because the range of sediment sizes that deposit in the delta model is between silt and medium sand.Using a shape factor of 40 results in K values that are within the range expected within a coastal delta (1 × 10 −4 -1 × 10 −8 m/s).In each model cell, the d 10 is determined using a linear interpolation with 1,000 values to calculate the cumulative density function representing the percent finer of the sediment deposited in that cell.Although other interpolation methods were considered (cubic, nearest), a linear interpolation was deemed most relevant by visual inspection of the percent finer graphs (Figure S2 in Supporting Information S1).K is calculated using a geometric mean of K ij contained within the shoreline.
Permeability (  ij ) for each cell in the delta plain is computed as a function of  ij , dynamic viscosity (   ), density of the pore water (   ), and gravitational constant (   ): where μ = 1.8 × 10 −3 kg/ms, ρ = 1,000 kg/m 3 , and g = 9.81 m/s 2 is used.The values for fluid density and viscosity assume fresh groundwater.Analysis of the salinity concentration in each model in the area of sediment deposition confirms this assumption.Permeability is calculated to provide an intrinsic property for the sediment alone, rather than K which depends on the fluid properties.The permeability (k) is calculated using a geometric mean of k ij .

Connectivity of High Permeability Bodies
A high permeability body is defined as any cluster of model cells that have k ij greater than the 75th percentile of the k ij within the model.The connectivity of the high permeability bodies (Co) is calculated by dividing the area of the largest high permeability body (Ahp) by the sum of the area of all high permeability bodies (Steel et al., 2022): The 75th percentile is used to define Co because it aligns with the qualitative analysis of the high permeability zones in each delta.We tested this method using the 90th percentile but found that the results were skewed by grid cells with anomalously high permeability-only a few grid cells in the tidal and fluvial deltas met this criterion.
Results from the 75th and 90th percentile were similar in the wave deltas.

Hydraulic Gradient
The slope of the delta plain, referred to herein as the topographic gradient of the delta (dE/dl), is calculated as an average topographic gradient between the initial channel and the shoreline.Elevation at the initial channel (Ec) is calculated as an average of the elevation in the 10 model pixels on both sides of the channel.This removes any impact that the topographical depression in the channel would have on the slope.A gradient is calculated for every pixel (n) along the shoreline by subtracting the elevation at the shoreline (Es n ) from Ec.An average slope is computed for each delta model by calculating the mean of the gradients: (5) The regional hydraulic gradient (dh/dl) is assumed to be equivalent to dE/dl.This is a reasonable assumption in a deltaic environment due to the gently sloping nature of the topography.Note, however, that local hydraulic gradients could be much greater in a real delta where groundwater pumping takes place.Furthermore, since we are only simulating the top Holocene sedimentation in modern deltas, we are interested in the shallow unconfined aquifer systems and not the deeper, confined, or semi-confined aquifer systems that exist as part of the transgressive and regressive nature of deltaic geologic history.

Specific Discharge
Specific discharge (or the groundwater flux; q) in each delta model is computed by multiplying the geometric mean of hydraulic conductivity (K) by the regional hydraulic gradient (dh/dl): Using the geometric mean of hydraulic conductivity (K) provides one representative q value for the entire delta and does not consider the complexity and spatial variability of the hydraulic conductivity within each delta.Additionally, we are only considering the horizontal groundwater flux through the delta since hydraulic conductivity only varies horizontally and is vertically uniform.A value of q computed in this manner also does not consider any recharge limitations that may exist and assumes that the system is topography-limited rather than recharge-limited.These simplifications mean that the computed values for q are likely not comparable to q in real deltas and therefore are only used as a tool to understand how differences in K and dh/dl may impact groundwater flow within the shallow subsurface in the three delta types.

Analysis of Sensitivity
We undertake three sensitivity analyses to determine how the morphodynamic influences, basin conditions, and geomorphic characteristics impact the hydrogeologic conditions that impact horizontal groundwater flow through the shallow subsurface of a delta.The sensitivity analyses are performed using a DGSA (Fenwick et al., 2014;Scheidt et al., 2018) to compute the sensitivity of the response to changes in the input parameters.The response parameters, or the hydrogeologic conditions of interest, include permeability (k), hydraulic gradient (dh/dl), and connectivity of the high permeability areas (Co).The response parameters were chosen because k and gradient 10.1029/2022WR034136 6 of 16 impact the groundwater flux.Co was selected because we were interested in the connectivity of the high permeability areas.The input parameters include fluvial discharge (Qav), wave height (Wa), tidal amplitude (Ti), sediment concentration (Cs), bathymetric gradient (Dgrd), median grain size (Dmm), the number of channels in the delta (Cn), delta shape (Sh), and shoreline rugosity (Sr).Models are sorted into a high and low class based on the response (Scheidt et al., 2018); for example, models that have a high k are sorted into a high response class while models that have a low k are sorted into a low response class when performing the sensitivity analysis for k.The sensitivity of the response to changes in an input parameter is quantified based on the distance between the prior frequency distribution and the frequency distribution computed for each class (Scheidt et al., 2018).A hypothesis test with a 95% confidence level (p < 0.05) is used to test the null hypothesis that the input parameter has no impact on the response (Scheidt et al., 2018).An influential, or critical parameter, rejects the null hypothesis while a non-influential, or insensitive parameter, fails to reject the null hypothesis.

Modeled Delta Geomorphology
Anderson et al. (2023) showed that the models used in this study can reproduce the geomorphologic characteristics expected in a fluvial, wave, and tidal delta (Figures 2a-2c).The modeled fluvial deltas have an average of seven channels, with only a few of these reaching from the delta apex to the shoreline (Figure 2a).The channels that do not connect to the delta apex are remnants of previous channelization that were deactivated, likely due to sediment backfilling or avulsion.The modeled tidal deltas have an average of eight channels, while the wave deltas have an average of two channels (Figure 2d).Although the fluvial and tidal deltas have a similar number of channels (Cn) when averaged across all delta models, tidal deltas have a wider range and overall increased channelization compared to the fluvial deltas.The distributary network in the tidal delta occurs throughout the delta plain, with many characteristic tidal channels that are wider at the shoreline and decrease in width inland (Figure 2c).In the wave delta, there is one main channel that cuts through the delta plain (Figure 2b).This main channel bifurcates near the channel outlet where a subaqueous mouth bar has formed (L.D. Wright & Coleman, 1973).Note that some of the channels in the fluvial and tidal deltas exhibit a blocky shape where the channel changes direction at a 90° angle; this is likely a result of channel avulsion within the model.Each of the three delta types can have as few as one channel (Figure 2d).For fluvial and tidal deltas, this happens more frequently in models that have a lower fluvial discharge (Qav), resulting in a smaller delta that experiences fewer avulsions and less tidal reworking.
Delta models primarily affected by fluvial influences have the smallest average shape (Sh = 1.1), indicating that the delta plain is elongated into the receiving basin (Figure 2e).Interestingly, the tidal deltas have the highest average Sh (1.8).Historically, wave deltas have been considered to have the highest Sh.However, due to the increased reworking of sediment from tides, many of the tidal deltas build sediment evenly across the southern boundary, resulting in a wider delta.In nature, many tidal deltas form in drowned valleys where the width of the delta is confined by the shape of the basin (J.Syvitski et al., 2022).This may explain why wave deltas are typically thought of as the widest deltas even if tidal reworking has the potential for lateral sediment transport.Additionally, the morphodynamic models only simulate waves that are perpendicular to the initial shoreline; this is often not the case in nature where wave incidence can be oblique to the river mouth, resulting in littoral drift that redirects sediment adjacent to the coastline (L.D. Wright, 1977).The average Sh for wave deltas is 1.5.Elongation of wave deltas into the receiving basin increases as fluvial discharge (Qav) increases.On average, wave deltas have the smoothest shoreline, ranging from 1.0 to 1.2.Tidal deltas have the most rugged shoreline as well as the widest range, with shoreline rugosity values ranging from 1.0 to 1.7.Some of the smaller values for Sr in tidal deltas can be explained by models that have some wave influence, reducing the overall tidal reworking and creating a smoother shoreline.

Permeability
The average permeability (k) within the delta models ranges from 3 × 10 −14 m 2 to 9.4 × 10 −11 m 2 (Figure 3a).The tidal deltas have the lowest k (1.3 × 10 −12 m 2 ), followed by the fluvial (2.7 × 10 −12 m 2 ) and the wave deltas (1.0 × 10 −11 m 2 ) (Figure 3a).The median k of the wave delta occurs near the top of the interquartile range, indicating that the distribution of k is negatively skewed.The percent of high, medium, and low permeability area in each delta type (Figures 3b-3e) indicates that, on average, wave deltas are comprised of 68% high permeability, 29% medium permeability, and 3% low permeability sediment.The fluvial deltas have 41% high and 33% low permeability material, with only 26% medium permeability.The tidal deltas have mostly medium permeability, comprising 53% of the delta material.Low permeability material makes up 33% of the tidal deltas and 14% is high permeability.
The spatial permeability maps (Figures 3f-3h) indicate that the higher permeability areas in the fluvial and tidal deltas occur in the middle of the delta and align with the distributary channel network.The higher permeability areas within the fluvial delta have k ≥ 2 × 10 −11 m 2 , with lower permeability material (k < 5 × 10 −13 m 2 ) along the shoreline.Most of the permeability within the tidal delta falls within 2 × 10 −13 m 2 to 7 × 10 −12 m 2 .The tidal delta has less low permeability material along the shoreline compared to the fluvial delta and has a few locations within the main channels where k ≥ 3 × 10 −10 m 2 .These higher permeability areas develop due to tidal influence of the river network.Modulation of the fluvial outflow results in lower water velocities and retention of course-grained material within the river network.
Less than 10% of the permeability in the wave delta is <1 × 10 −12 m 2 (Figure 3g).The regions with the greatest permeability in the wave delta (k ≥ 2 × 10 −11 m 2 ) occur along the sides of the main channel, where the delta has built channel margins.Within the channel margins, some linear high permeability structures appear; it is possible that these are swash bar features, indicative of wave influence.The part of the delta plain located adjacent to the channel margins is characterized by a smooth shoreline and medium permeability (k < 2 × 10 −11 m 2 ).This part of the delta plain exhibits little variation in k because most of the sediment that deposits adjacent to the channel margins is fine sand.Smaller grain sizes are unable to deposit within the delta plain because increased wave action leads to higher water velocities, inhibiting the deposition of fine-grained material.However, the wave action in these models is not powerful enough to carry the coarse-grained sediment in bedload transport away from the river mouth.

Connectivity of High Permeability Bodies
The maximum number of connective bodies in each of the delta types is 76 for fluvial deltas, 27 for wave deltas, and 62 for tidal deltas.The median connectivity (Co) in all the delta types is 0.99, indicating that most of the high permeability bodies are connected (Figure 4a).Wave deltas have the largest interquartile range for connectivity (Co = 0.75 to 1.0), implying that the connective bodies within wave deltas can be more discontinuous than in fluvial and tidal deltas.This results from the connective bodies that exist within the channel margins being separated by lower permeability material within the channel.In wave delta models that have large fluvial discharge (Qav), the high water velocity in the channel is too great for fine-grained sediments to deposit, resulting in higher permeability within the channel and higher connectivity between the channel margins.When Qav is small compared to wave height (Wa), lower water velocities reduce the transport competence of the river and allow for finer sediments to settle within the channel, separating the high permeability bodies.The increased variability of sediment deposition within the channel compared to the channel margins appears to reduce the overall permeability within the channel (Figure S3 in Supporting Information S1).In general, the high permeability bodies within the wave delta are larger than those in the fluvial or tidal deltas and are contained within the channel margins rather than spread throughout the delta.The median size of a connective body in a wave delta is 74% of the total delta area, compared to 32% of the total delta area in both the fluvial and tidal deltas.
The percent of the channel network that contains a high permeability body is 45% in fluvial deltas, 35% in tidal, and 22% in wave (Figure 4c).The high permeability bodies within the fluvial delta correspond to the active distributary channels as well as in previously channelized areas, located in the middle of the delta plain (Figure 4d).The high permeability bodies in the wave delta are less contained within the channel network than in fluvial and tidal deltas because the high permeability bodies exist adjacent to the channel where channel margins are constructed (Figure 4e).Although most of the high permeability bodies in the tidal delta align with the tidal network, some of the tidal channels that exist outside of the active delta (channels that form near the southern boundary) do not contain high permeability material (Figure 4f).

Hydrogeologic Properties
The median hydraulic gradient (dh/dl) in the delta models is 2.9 × 10 −4 m/m (Figure 5a).The dh/dl is largest in the tidal delta, with a median gradient of 1.4 × 10 −3 m/m.The increased dh/dl in the tidal delta results from an increase in accommodation space during flood tides, allowing sediment deposition at higher elevations.The median dh/dl in the fluvial and wave deltas are 2.9 × 10 −4 m/m and 1.4 × 10 −4 m/m, respectively.
Although there is variation in hydraulic gradient (dh/dl) between the different delta types, there is little variation between the medians of specific discharge (q; Figure 5b).Wave deltas have the greatest median q (1.4 × 10 −8 m/s), while tidal deltas have the lowest q (7.2 × 10 −9 m/s).Although tidal deltas have the largest dh/dl, the lower average hydraulic conductivity (K) results in a smaller q compared to the wave and fluvial deltas.Similarly, although the wave delta has the smallest dh/dl, the increased K results in a larger q in the wave delta models.

Sensitivity Analysis
Permeability is most sensitive to changes in the tidal amplitude (Ti) and wave height (Wa), as well as sediment concentration (Cs) and the number of channels (Cn) (Figure 6a).The sensitivity of permeability (k) to changes in Ti and Wa is explained by the ability of the delta to retain coarse-grained material within the delta.Deltas with high Wa prohibit the deposition of fine-grained material and retain coarse-grained material within the channel margins.In comparison, larger Ti results in more sediment of medium permeability that deposits around the channel network, which is extensive throughout the delta plain.The distinction between morphodynamic influences within wave and tidal deltas is reflected in the number of channels within a delta, making k also sensitive to changes in Cn.The sensitivity of k to Cs indicates that the grain size that is deposited in the delta is dependent on the amount of sediment coming into the basin.However, k is insensitive to changes in the median grain size (Dmm).
The hydraulic gradient (dh/dl) within a delta is also sensitive to changes in the tidal amplitude (Ti), in addition to the number of channels (Cn) and shoreline rugosity (Sr) (Figure 6b).Gradient is highly sensitive to Ti because tides can deposit sediment at higher elevations when the water level rises.Once again, an increased Cn is indica tive of tidal influence, resulting in dh/dl being sensitive to Cn.The dh/dl is also sensitive to changes in Sr, which is better able to distinguish fluvial deltas from tidal deltas.
Connectivity of high permeability material (Co) is most sensitive to changes in fluvial discharge (Qav), wave height (Wa), and tidal amplitude (Ti); it is not sensitive to any of the geomorphic characteristics (Figure 6c).The sensitivity of Co to Qav and Wa may be explained by the relative proportion of fluvial and wave power.When fluvial and wave power are similar, the convergence of fluvial discharge and wave action results in lower water velocity, allowing fine sediment to settle and create a barrier between the two high permeability channel margins.When either fluvial or wave power dominates within the channel, the water velocity remains high and prohibits the deposition of fine sediment within the channel.

Discussion
Many previous groundwater modeling studies focusing on coastal deltas have used 1D or 2D models that evaluate vertical groundwater flow through a column or along a longitudinal cross-section of the delta (Bridger & Allen, 2006;Delsman et al., 2014;Larsen et al., 2017;Van Pham et al., 2019).The lack of spatially varying lateral heterogeneity and anisotropy within deltas results in a limited understanding of how permeability structures vary across the delta plain and how this impacts horizontal groundwater flow through a delta.We focus on quantifying and characterizing the horizontal spatial variability in the deltaic subsurface that arises because of varying morphodynamic influences.Our findings that deltaic permeability and hydraulic gradient are sensitive to changes in morphodynamic influences, as well as variations in the channel network and shoreline rugosity, are promising because they suggest that some knowledge of the subsurface properties of a delta can be inferred by simply identifying if the delta is fluvial, wave, or tidal.The quantification of permeability, hydraulic gradient, connectivity of permeable bodies, and the location of high permeability bodies for three end-member delta types provides a general understanding of what the hydrogeologic system in a delta might look like spatially without having to conduct costly field campaigns.Our results also support previous work that has linked permeable bodies and preferential pathways to current and previous channelization (Kolker et al., 2013;Steel et al., 2022;Xu et al., 2021).
The permeability maps generated herein are limited by the cell size within the modeling domain, the steady-state boundary conditions, and homogeneity in the vertical direction.To reduce runtimes of the morphodynamic models, a cell size of 200 m × 200 m is used; this approach assumes that each grid cell is homogeneous and isotropic, neglecting horizontal heterogeneity that varies within 200 m.Since deltas exhibit variability within a 200 m spatial scale, the large cell size is a limitation of this work.Overbank deposits and high permeability features caused by small channels (less than 200 m wide) are two features that are not reflected in the permeability maps due to the model setup.Of these two features, overbank deposits are more likely to have an impact on groundwater flow through a delta since these deposits typically have low permeability and may act as a barrier to horizontal groundwater exchange between the permeable channels/paleochannels and the surrounding sediment.Using a smaller cell size may allow for overbank deposits to be more distinguishable within the permeability maps by adding more low permeability material around the channel networks.The formation of overbank deposits requires flooding that is often associated with extreme weather events or seasonal high flows; using steady-state boundary conditions means that the variability required to form overbank deposits was not taken into account.Although overbank deposits are features in many deltas, Michael and Voss (2009) noted the absence of overbank deposits in the Bengal Basin and attributed this to the mobile nature of the channels in the tidal environment.The results of this study would benefit from the inclusion of seasonal variability within the boundary conditions and from a smaller spatial resolution that can account for variability at smaller spatial scales.
The vertical structure of sediment deposition is not recorded in these models.This means that any sediment contained within a cell is assumed to be "well-mixed," producing vertically homogenous sediment deposits within the model.This limitation means that a vertical sediment structure cannot be discerned.Channel stacking patterns are important when considering the vertical connectivity of fluvial sediment deposition within deltas around the world.Further, the ability for channels to erode through marine confining layers and create vertical discontinuities is critical to account for when evaluating regional-scale subsurface structure in a delta.Improvements to this work may include the use of 3D modeling or implementation of the stratigraphy module in Delft3D (Deltares, 2013), allowing for the preservation of sedimentation records vertically throughout the model simulation.Using geometric averaging, we provide a first approximation of the rate at which groundwater moves through deltas around the world.

Permeability and Morphodynamics
We estimate that overall permeability in deltaic landforms has a median value of 4.0 × 10 −12 m 2 , relating to a hydraulic conductivity value of 2.14 × 10 −5 m/s.Across all delta models, 42% of permeabilities represent medium sand (k > 2.7 × 10 −10 m 2 , K > 5 × 10 −5 m/s), 40% fine sand (2.7 × 10 −10 m 2 ≥ k > 1.1 × 10 −11 m 2 , 5 × 10 −5 m/s ≥ K > 2 × 10 −6 m/s), and 18% silt (1.1 × 10 −11 m 2 ≥ k > 5.5 × 10 −15 m 2 , 2 × 10 −5 m/s ≥ K > 1 × 10 −9 m/s).These values are consistent with the permeability values that deltaic groundwater flow and solute transport modeling studies have used (Table 1).Additionally, the simulated thickness of the sediment within the delta models ranged from 15 to 25 m near the coast (Figure S4 in Supporting Information S1), which is comparable to the reported thickness of shallow aquifers in some large coastal deltas including the Vistula and Yangtze (Anderson, 2023;Cao et al., 2013;Jaworska-Szulc, 2009;Li et al., 2006).The percent of sand and fine sediment in the Mississippi Delta reported by (Williamson et al., 1990) closely resembles the proportion of high, medium, and low conductivity areas for fluvial deltas herein (Table 1, Row 1).Most deltas that have in-depth groundwater studies are tidal (Table 1), which may reflect an overall lower permeability, considering our finding that tidal deltas have the lowest average permeability and the greatest proportion of medium and low permeability compared to the overall delta size.
Many groundwater modeling studies report the existence of clays and peat layers (Larsen et al., 2017;Michael & Voss, 2009;Minderhoud et al., 2015;van Engelen et al., 2019;Van Pham et al., 2019).However, lower permeabilities that would be representative of clays (k < 5.5 × 10 −15 m 2 , K < 1 × 10 −9 m/s) were not present in any of our delta models.Analysis of sediment profiles within the delta (Figure S3 in Supporting Information S1) indicates that clays and some silts do not deposit within the delta plain due to increased water velocities.Clays and silts are found in the subsurface of deltas around the world and are often derived from marine transgression and seasonal inundation.Clay and silt units deposited during periods of sea level rise often make up the aquitard structures in coastal deltas and can be barriers to vertical groundwater flow.The results of this study would benefit from including temporal variation in the boundary conditions and simulating seasonal flooding events that lead to the deposition of more fine-grained material within the delta.The addition of sea level change within the modeling structure and the inclusion of vertical variability in the sediment deposits would likely increase the preservation of fine-grained sediments within the subsurface and would allow for analysis of subsurface permeability vertically.However, we do find that our overall permeability is within the range of permeability found in natural deltas.One explanation for the absence of clay within our models may be that the models are forced with constant boundary conditions.Overbank deposits require flooding, often associated with extreme weather events or seasonal high flows, were not accounted for in this modeling.More accurate spatial representations of deltaic permeability may be acquired by using temporally varying forcings to simulate the formation of overbank deposits and their preservation in end-member morphodynamic environments.
Permeability (k) representative of coarse-grained sands (k > 1.1 × 10 −9 m 2 , K > 6 × 10 −3 m/s) is also not present in the delta modeling even though grain sizes up to 2 mm are included in the incoming sediment within the models.Analysis of the sediment distributions indicates that the coarsest-grained material does deposit within the river channels in the tidal and fluvial deltas, and in the channel margins in the wave deltas (Figure S3 in  Supporting Information S1).However, the proportion of finer-grained sands that also deposit within these model cells is greater, resulting in a hydraulic conductivity that is more representative of the fine-grained sand material.
The permeability averages reported in this study are not sensitive to changes in the median incoming sediment grain size within the model, suggesting that deltaic permeability is not controlled by the grain size of sediment being delivered to the delta.This is especially true in wave deltas where wave action increases the water velocity and inhibits the deposition of fine-grained sediments within the delta plain and carries the small particle sizes offshore (Wright & Coleman, 1973).Much of the fine-grained sediments found in real deltas around the world are a product of sea level variation that changes a near-shore environment into an offshore environment (Tanabe et al., 2006).Here we only consider the sediment deposition that occurs within the deltaic basin during delta formation (when sea level regression slows).Therefore, permeability within the delta is dictated by the presence of morphodynamic influences.This likely results from the natural sorting abilities of the morphodynamic influences within the delta; high water velocities arising from increased fluvial discharge, waves, and tides prohibit the deposition of small grain sizes (Wright & Coleman, 1973).However, these influences promote bedload transport of the larger grain size material.Lower water velocities, specifically outside the impact of the fluvial discharge and in the absence of waves and tides, allow for more deposition of fine-grained sediments, which results in lower permeability.The resulting permeability distributions appear differently in each delta type based on which influences are present within the model and the strength of those influences throughout the receiving basin.

Groundwater Flow
We found that despite the variation in the average permeability in the three delta types, specific discharge remains relatively similar across the deltas due to variations in the hydraulic gradient.The increased permeability in the wave deltas is offset by a smaller hydraulic gradient due to low topographic relief.On the other hand, the increased hydraulic gradient in the tidal deltas stemming from a larger topographic relief is counteracted by overall lower permeability.Although variations in the specific discharge of the three deltas types is small, delta type is important when considering barriers to groundwater flow.Steel et al. (2022) state that the presence of low permeability barriers may be more important to groundwater flow through the delta than the permeable pathways; this conclusion was made by physically modeling the formation of a delta within a basin that is void of waves and tides.The permeability maps in this study confirm that fluvial deltas have an accumulation of low permeability material along the shoreline that blocks some of the permeable pathways from extending from the delta apex to the shoreline.This may also be relevant for wave and tidal deltas that have an accumulation of fine-grained sediment adjacent to the delta plain or offshore in the subaqueous delta.
Our models show a close link between the high permeability bodies within a delta and the location of the channel network, which is consistent with other generic delta modeling studies that have estimated permeability or hydraulic conductivity changes laterally (Hariharan et al., 2021;Steel et al., 2022;Xu et al., 2021).The high permeability material within fluvial and tidal deltas is associated with current and paleochannelization and is primarily located within the channels.Furthermore, the high permeability bodies are highly connective and extensive throughout the delta.The high connectivity of permeable bodies horizontally was also found by Xu et al. (2021) and Steel et al. (2022) with the numerical and physical modeling of a fluvial delta.Modeling the formation of deltas using numerical morphodynamic models can preserve the high permeability features that are created due to previous channelization which is not visually apparent on the delta surface due to avulsions.This is particularly important when considering permeability structures in fluvial deltas, as more of the high permeability corresponds to channels that have been backfilled.
The high permeability features in the wave delta are not directly located within the channels but rather on the channel margins.Many wave delta models exhibit lower permeability within the channel network.Although there is more coarse sand deposited within the channel compared to the channel margins, there is also increased fine sand.The increased variability of sediment deposition within the channel compared to the channel margins appears to reduce the overall permeability within the channel.The blob-like high permeability bodies located within the channel margins in the wave delta are less likely to be connected compared to the high permeability bodies that follow the channel network in the fluvial and tidal deltas.Although detailed spatially varying permeability maps of many real-world wave deltas are often not available, similarities can be drawn between the modeled and real wave deltas.High permeability swash bar features are known to form in wave deltas and are apparent in the Sao Francisco and Senegal deltas (Wright, 1977) and may explain the high permeability linear features contained within the channel margins.The balance of the morphodynamic influences shaping deltas has been subject to change in the recent past and will continue to change in the future.The effects of this are already evident, with river damming restricting the amount of discharge and sediment entering deltas globally (Giosan et al., 2014;Nienhuis et al., 2020).As sea level rises, it is also likely that wave and tidal action will increase, shifting the balance toward marine dominance in many global deltas (Nienhuis et al., 2020).Increased extreme weather events combined with rising sea levels not only threaten to inundate deltaic lands but will likely result in the deposition of more fine-grained surface material intermixed with highly erosive episodes.The shallower elevation and hydraulic gradient in wave deltas make these landforms more susceptible to inundation as well as groundwater salinization through marine encroachment.Additionally, the high connectivity of the deltaic subsurface may allow for faster salinization of the subsurface through permeable pathways.

Conclusions
Quantitative evaluation of groundwater supply in coastal deltas requires the estimation of large-scale hydrogeologic characteristics that control groundwater flow.This study evaluated the horizontal permeability in fluvial, wave, and tidal deltas.The sensitivity analysis suggests that permeability, hydraulic gradient, and the connectiveness of high permeable bodies can be characterized through the identification of morphodynamic influences (fluvial, wave, and tidal) and geomorphic characteristics (number of channels and shoreline rugosity) within a delta.Using numerical morphodynamic modeling of delta formation in 171 unique models, we found that: 1.The overall permeability in deltaic landforms has a median value of 4 × 10 −12 m 2 , relating to a hydraulic conductivity value of 2.14 × 10 −5 m/s.The average hydraulic gradient is 3.9 × 10 −4 .Wave deltas are the most permeable but have the smallest hydraulic gradient while tidal are the least permeable and have the highest hydraulic gradient.2. The high permeability bodies are associated with current and previous channelization in a delta and are highly connected horizontally in the shallow subsurface.This high connectivity potentially allows for salinization through permeable pathways as sea level rises.The channel network is most evenly distributed in tidal deltas, resulting in high permeability bodies located throughout the entire delta plain.High permeability bodies in wave deltas are only located in the channel margins and are not pervasive throughout the entire delta.3. Wave deltas may be most susceptible to inundation and groundwater salinization through marine encroachment due to the smaller hydraulic gradient and increased permeability.
by first finding the K of each individual cell (200 m × 200 m) within the model domain,  = Cd10 2 .

Figure 2 .
Figure 2. Examples of a modeled (a) fluvial, (b) wave, and (c) tidal delta.The (d) number of channels in the delta (Cn), (e) delta shape (Sh), and (f) shoreline rugosity (Sr) for the three delta types and all delta models.

Figure 3 .
Figure 3. (a) The geometric mean of the permeability (k) across the entire delta (all model cells) and (b-e) the portion of high, medium, and low permeability material in fluvial, wave, tidal, and all delta models.Spatial permeability map of a (f) fluvial, (g) wave, and (h) tidal model.

Figure 4 .
Figure 4. (a) The connectivity of high permeability bodies (Co), (b) size of the largest permeable body normalized by delta area, and (c) percent of the channel network that contains a high permeability body in fluvial, wave, tidal, and all delta models.Maps of the high permeability bodies for a (d) fluvial, (e) wave, and (f) tidal model.
The delta type reflects the dominant morphodynamic influence controlling the geomorphic characteristics in the delta (based on the original data published inNienhuis et al., 2020).a Values from the Holocene layer only.

Table 1
Horizontal Hydraulic Conductivity Values Used in Deltaic Groundwater Modeling Studies