Characterizing Offshore Freshened Groundwater Salinity Patterns Using Trans‐Dimensional Bayesian Inversion of Controlled Source Electromagnetic Data: A Case Study From the Canterbury Bight, New Zealand

The study of offshore freshened groundwater (OFG) is gaining importance due to population growth and environmental pressure on coastal water resources. Marine controlled source electromagnetic (CSEM) methods can effectively map the spatial extent of OFG systems using electrical resistivity as a proxy. Integrating these resistivity models with sub‐surface properties, such as host‐rock porosity, allows for estimates of pore‐water salinity. However, evaluating the uncertainty in pore‐water salinity using resistivity models obtained from deterministic inversion approaches presents challenges, as they provide only one best‐fit model, with no associated estimate of uncertainty. To address this limitation, we employ trans‐dimensional Markov‐Chain Monte‐Carlo inversion on marine time‐domain CSEM data, acquired in the Canterbury Bight, New Zealand. We integrate resistivity posterior probability distributions with borehole and seismic reflection data to estimate pore‐water salinity with corresponding uncertainty estimates. The results highlight a low‐salinity groundwater body in the center of the survey area, hosted by consecutive silty‐ and fine‐sand layers approximately 20–60 km from the coast. The posterior probability distribution of resistivity models indicates freshening of the OFG body toward the shoreline within a permeable, coarse‐sand layer 40–150 m beneath the seafloor, suggesting an active connection between the OFG body and the terrestrial groundwater system. The approach demonstrates how Bayesian inversion constrains the uncertainties in resistivity models and subsequently in pore‐water salinity estimates. Our findings highlight the potential of Bayesian inversion to enhance our understanding of OFG systems and provide uncertainty constraints for hydrogeological modeling, thereby contributing to sustainable water resource development.


Introduction
The rising demand for potable water in coastal regions has promoted growing research interest in detecting offshore freshened groundwater (OFG) worldwide (Post et al., 2013).OFG emplacement is commonly explained through either meteoric recharge during sea-level low-stands (e.g., Cohen et al., 2010;Meisler et al., 1984;Person et al., 2003) and/or through active meteoric recharge via permeable connections between onshore and offshore aquifers (e.g., Hong et al., 2019;Johnston, 1983;Weymer et al., 2022).Due to the ongoing pressure on terrestrial aquifers caused by environmental and anthropogenic stressors, offshore groundwater research has gained new urgency.One proposition is that near-coastal OFG reservoirs can serve as unconventional sources of drinking water, especially along densely populated coastal regions and on islands, where terrestrial aquifers cannot maintain a sufficient supply for inhabitants, agriculture, and industry (Bakken et al., 2012;Micallef et al., 2021;Zamrsky et al., 2022).
Evidence of OFG has been documented in coastal embayments and continental shelves worldwide (e.g., Attias et al., 2020;Cambareri & Eichner, 1998;Gustafson et al., 2019;Haroon et al., 2021;Haroon, Lippert, et al., 2018;Hathaway et al., 1979;King et al., 2022;Martin et al., 2007;Micallef et al., 2020;Thomas et al., 2023), demonstrating that low-salinity groundwater within the seafloor (water depths <1,000 m) occurs on a global scale (Post et al., 2013).Global OFG volume has been estimated at 10 6 km 3 (Micallef et al., 2021), yet these estimates are poorly constrained since they are based on first-order assumptions derived from point-source borehole measurements and hydrogeological modeling.A general lack of integrated geophysical data and hydrogeological investigations along continental shelves precludes a seamless accounting for more precise volume and salinity estimates on both global and regional scales (Weymer et al., 2022).In order to argue for the sustainable use of OFG as an unconventional resource, volumes, geometries, salinities, as well as environmental impacts need to be precisely quantified before exploitation.
The exploration of OFG requires a non-invasive and non-destructive geophysical mapping tool that is sensitive to pore-water salinities.Since the only physical material property sensitive to pore-water salinity is electrical resistivity, marine controlled-source electromagnetic (CSEM) studies have recently gained attention in the realm of OFG research.Several applications of marine CSEM methods have proven effective in identifying OFG residing in various offshore geological settings (e.g., Attias et al., 2020;Blatter et al., 2019;Evans, 2007;Gustafson et al., 2019;Haroon et al., 2021;Haroon, Lippert, et al., 2018;King et al., 2022;Levi et al., 2018;Lippert & Tezkan, 2020;Micallef et al., 2020).The electrical resistivity distribution of the seafloor sediments, which can be obtained from CSEM data, is controlled by the electrical conductivity (1/resistivity) of the pore-fluid, the porosity of the host rock, the degree of its cementation, and the surface electrical conductivity (Archie, 1942;Hoefel & Evans, 2001).Thereby, resistivity models are utilized as a proxy to delineate the spatial extent of OFG systems along continental shelves by applying the empirical Archie's relationship (Archie, 1942) which relates the resistivity of the seafloor to the resistivity of the pore-fluid (e.g., Gustafson et al., 2019;King et al., 2022;Micallef et al., 2020).Currently, only one study assesses the uncertainty of the derived resistivity models and their impact on the uncertainty of OFG reservoir characterization (Blatter et al., 2019).This assessment is particularly important when aiming to validate hydrogeological models with geophysical data, an approach that has been identified as instrumental in the future assessment of OFG resource potential, evolution, and responses to external drivers (Arévalo-Martínez et al., 2023;Weymer et al., 2022).
From a geophysical perspective, one of the main challenges in the assessment of OFG bodies is to reliably predict the pore-fluid salinity distribution based on electrical resistivity models.Two primary factors contribute to uncertainty in the OFG salinity estimation: (a) the uncertainty of the best-fit resistivity model and (b) unknown variations in lithology-dependent parameters in Archie's relationship.The majority of CSEM studies apply a smoothness constraint in a deterministic inversion approach and subsequently estimate pore-fluid salinities by assuming representative values for the cementation factor and host-rock porosity in Archie's relationship (Archie, 1942).In most cases, the comparison between predicted salinity distributions derived from geophysical models and hydrogeological modeling or ground-truth measurements at boreholes achieves only a qualitative match (e.g., Gustafson et al., 2019;Micallef et al., 2020), possibly because model parameter uncertainties are not accounted for in the analysis.A more quantitative assessment of OFG salinity distributions can be achieved if uncertainties are addressed within the entire evaluation chain of deriving pore-water salinity models from integrated CSEM analysis.We adopt an approach similar to that of Blatter et al. (2019) for integrating uncertainties in porosity values and Archie's parameters when estimating pore-water salinity.A notable distinction in our approach is that instead of relying solely on in situ porosity measurements obtained from boreholes to construct porosity distributions, we incorporate porosity-depth profiles extracted from interpreted seismic sections.This approach proves to be practical, particularly when offshore borehole measurements are scarce or entirely unavailable.
An uncertainty evaluation of pore-water salinity predictions derived from CSEM measurements will allow inference about the genesis of OFG and its potential connection to their terrestrial baseline, as well as provide a cornerstone for future monitoring of freshwater/saltwater interfaces with electromagnetic methods.These uncertainty estimates ideally require a statistical analysis of all parameters in Archie's empirical relationship (Archie, 1942), that is, bulk resistivity, porosity, and cementation factor.Deterministic CSEM inversion approaches use an objective function that is constrained by additional regularization terms such as a roughness constraint (e.g., Constable et al., 1987) to iteratively optimize the data misfit between measured and modeled data and ensure numerical stability.These deterministic inversion approaches are computationally efficient when searching for one best-fit model (e.g., Abubakar et al., 2008;Constable et al., 1987) and are typically used for 2-D and 3-D inversion of CSEM data.Yet, they cannot quantify the inherent uncertainty in model parameters.Moreover, smoothed resistivity models pose additional challenges when integrated with other geophysical data, such as multi-channel seismic (MCS) reflection or borehole data, due to disparate spatial resolution.However, integrating electrical resistivity models with lithological information derived from seismic/acoustic methods is essential for the regional mapping of OFG (e.g., Weymer et al., 2022), as lithological sequences and geological heterogeneities govern the distribution of OFG (Zamrsky et al., 2020).
In contrast to deterministic inversion approaches, Bayesian inversion algorithms treat parameters as random variables that are probabilistically sampled, resulting in a posterior probability distribution of resistivity models that fit the data, instead of a single "best-fit" model.These distributions can be used to provide rigorous estimates of model parameter uncertainties (Malinverno, 2002); however, the large computation time associated with Bayesian inversion currently limits their application to 1-D modeling.We apply a trans-dimensional Bayesian inversion algorithm to our CSEM data, which provides layering structure estimates and includes the number of model parameters as unknown in the inversion.This inversion approach has been applied to various geophysical data (e.g., Bodin & Sambridge, 2009;Dettmer et al., 2011;Gallagher, 2012;Gehrmann et al., 2015;Ray et al., 2013).
The goal of this study is to introduce a workflow using a 1-D Bayesian inversion to evaluate uncertainty in pore-water salinity predictions derived from integrated analysis based on CSEM, MCS, and borehole data.We base our analysis on a data set acquired in the Canterbury Bight, NZ (cf., Figure 1), where all required data sets are available and an OFG body is known to exist within the shallow seafloor (Micallef et al., 2020).First, we validate the developed workflow at the site of borehole U1353 drilled International Ocean Discovery Program (IODP) Expedition 317 (Fulthorpe et al., 2011; Figure 1) for ground-truthing and subsequently apply it at locations further away from the borehole, where the aquifer response detected by CSEM is greatest.Lastly, we discuss insights gained from this analysis regarding the evolution of Canterbury Bight's OFG system.

Regional Setting and Available Data
Canterbury Bight is located along the eastern coast of the South Island of New Zealand.The continental margin is constrained on its western side by the Southern Alps (>3,500 m height).The continental shelf, spanning an approximate area of 13,000 km 2 , consists of a ∼1 km-thick progradational successions of shelf-slope deposits (e.g., Browne & Naish, 2003;Lu & Fulthorpe, 2003).The Canterbury plains are underlain by a Quaternary succession of cyclically stacked fluvio-deltaic gravel, sand, and mud (Browne & Naish, 2003).These sediments were eroded from the Southern Alps during glacial periods, carried by high energy braided rivers, and deposited eastwards onto the subsiding (0.2-0.5 m ka 1 ) Canterbury margin (Browne & Naish, 2003;Rowan et al., 2012).To determine the relative contributions of global sea level fluctuations (eustasy) and local tectonic, sedimentary and oceanographic processes in controlling continental margin depositional cyclicity, three IODP sites were drilled along a coast perpendicular transect during IODP Expedition 317 (Fulthorpe et al., 2011;Figure 1a).A subtle decrease in pore-water salinity was observed at Site U1353, which is located closest to the shoreline at a water depth of 85 m and a distance of about 45 km from the coast.This salinity anomaly was identified at a depth between 20 and 60 m below the seafloor (bsf) and contains brackish pore-water with a salinity of approximately 24 practical salinity units (psu), in contrast to the background seawater salinity of 34 psu in the survey area (Fulthorpe et al., 2011).
In 2018, a joint marine CSEM and seismic survey was conducted during expedition TAN1703 (Mountjoy et al., 2017) within the framework of the MARCAN project (Micallef et al., 2018).The survey covered four profiles spanning a total length of 175 km across the Canterbury Bight and utilized a time-domain, seafloor-towed CSEM system.This system employs a 100 m long horizontal electric transmitter dipole (Tx) which is powered by a current-controlled signal generator hosted inside a decompressor ("pig") towed in front of the seafloor array.The signal generator has the capability to transmit currents up to 20 A. The transmitted source signal has a square waveform with a 50% duty cycle and a period of 4 s.Four inline electric receiving dipoles are connected to autonomous data loggers (Rx1-Rx4) which are towed behind the transmitter dipole at offsets of 150 m, 250 m, 400 m, and 650 m measuring the inline electric field components at 10 kHz (Gehrmann et al., 2016).The pig also contains a conductivity, temperature, and depth (CTD) probe to measure water depth and seawater resistivity, and an acoustic transponder to allocate the position of the seafloor system.Time-domain CSEM data were obtained at 267 stationary waypoints (WP) spaced every 500 m along lines 2, 4, and 5, and every 1 km along line 7. MCS data was acquired using a GeoEel Digital seismic streamer (Geometrics).The streamer consisted of 24 channels with a group interval of 12.5 m, and offsets between 22 and 308 m.The seismic source was a single mini GI-gun (13/35 cubic inch) towed at 1.5 m below the sea-surface.
Resistivity models derived from a 2-D deterministic inversion of the acquired CSEM data are presented by Micallef et al. (2020) using a time-domain extension of MARE2DEM (Haroon, Hölz, et al., 2018;Key, 2016).Moreover, Micallef et al. (2020) interpreted the MCS reflection profiles in terms of geological facies that include clay, silt, fine sand, coarse sand, and gravel, respectively.For conducting hydrogeological modeling, Micallef et al. (2020) assigned porosities of 45% to clay and silty material, 40% to fine sands, 35% to coarse sands, and 20% to gravel.However, it is important to note that these values represent first-order estimates, and their uncertainty is likely to increase at greater distances from the borehole U1353.This is due to limitations in the interpretation bias in seismic facies classification, inadequate information on sediment sorting and compaction, and local heterogeneities.
According to Micallef et al. (2020), the most prominent resistivity anomaly (>25 Ωm) is identified north of the borehole on the coast perpendicular line 4 and along the coast parallel line 7 (cf., Figure 4 in Micallef et al., 2020) This anomaly is significant as marine sediments in the region typically exhibit resistivity values ranging from 0.3 Ωm to 3 Ωm when fully saturated with seawater (Fulthorpe et al., 2011).Therefore, the identified resistivity anomalies suggest the presence of a large OFG body to the north of the borehole.The results of Micallef et al. (2020) also indicate that site U1353 serves as the southern boundary of a substantial OFG system that extends northwards up to line 5 (Figure 1).Micallef et al. (2020) estimated pore-water salinities from the 2-D resistivity models by applying Archie's relationship assuming constant porosities of 20, 30, and 40% to define plausible upper and lower limits of the potential OFG volumes.Their results suggest that the volume of OFG may increase by a factor of four as porosity increases from 20 to 40%.The applicability of Archie's equation in the survey area is supported by the interpretation of seismic sections, which reveals that the subsurface is primarily composed of sand, and clay content confined to thin surface layers (Micallef et al., 2020).
To enhance the accuracy of the salinity estimates and address associated uncertainties, we selected waypoints on line 2, intersecting with the borehole, and on line 4, which features the identified resistivity anomaly, for a Bayesian analysis.The workflow of this analysis is explained in the following section.

Methods
Marine time-domain CSEM surveys typically use a horizontal electric transmitter dipole to induce electrical currents into the surrounding seawater and seafloor.The electromagnetic (EM) signals propagate through the conductive seawater, where they are quickly attenuated, and through the more resistive seafloor.The Earth's responses are recorded by one or more receivers at different offsets from the transmitter.
To quantify model parameters with corresponding uncertainty estimates using CSEM data, we employ a 1-D trans-dimensional Bayesian inversion using a Markov Chain Monte Carlo (MCMC) sampling approach (e.g., Mosegaard & Tarantola, 1995;Sambridge & Mosegaard, 2002) and the Metropolis-Hastings-Green acceptance criterion (Green, 1995;Hastings, 1970;Metropolis et al., 1953).The term "trans-dimensional" denotes that transitions between models of different dimensions are enabled, that is, the number of sub-seafloor layer interfaces, and therefore the number of associated model parameters (resistivity and layer thickness) is variable.Here, we are adopting the algorithm from Gehrmann et al. (2015) that has previously been developed to address the geo-acoustic problem (Dettmer et al., 2010).Following the concepts of previous publications (Blatter et al., 2019;Gehrmann et al., 2015;Minsley, 2011;Ray & Key, 2012), we compute the posterior probability density (PPD) of resistivity models that are constrained by measured CSEM data and independent prior information about the site.These probabilities are then converted to probabilities of pore-water salinity by applying Archie's relationship (Archie, 1942).

Bayesian Formulation
Here, our Bayesian inversion employs MCMC sampling to estimate the PPD.The computational cost for Bayesian inversion can be large and therefore limits its applications to 1-D layered resistivity models that are computationally less expensive compared to 2-D/3-D modeling.However, the predominantly horizontal stratigraphy observed in the CSEM experiment's scale of our working area allows for a 1-D inversion.
The PPD combines prior knowledge about the model parameters and data information (e.g., Gelman et al., 2014), as expressed by Bayes' theorem, where d and m are the vectors of data and model parameters, respectively, and p represents the probability density functions (PDF).The PPD is the conditional probability of model parameters given the observed data, p(m|d).
The term p(d|m) is the conditional probability of the data given the model, which is interpreted as the likelihood of the model, L(m), for observed data.The prior p(m) is a probability density function of model parameters and is independent of the data.In this study, uniform distributions with parameter bounds are chosen as priors so that the solution is primarily constrained by the data.Specifically, we constrain the layer thickness and bulk resistivity values to fall within the respective ranges of 2-700 m and 0.2 to 200 Ωm.The resistivity limits are based on the resistivity range observed in the 2-D models.Additionally, this bulk resistivity range is compatible with a porefluid salinity range of 0.1 psu to values greater than 34 psu, given the expected porosity range in the region.The probability p(d) is the Bayesian evidence and normalizes the probability to unity.It can be ignored in this work (e.g., Sambridge et al., 2006).
For the likelihood function L(m), we assume that data errors are zero-mean and Gaussian distributed: where data misfit χ 2 quantifies the fit of predicted data f(m) and observed data d with covariance matrix C d .

Water Resources Research
10.1029/2023WR035714 (3) In this study, we use 1-D trans-dimensional MCMC algorithm with parallel tempering (Dettmer & Dosso, 2012) to efficiently sample the parameter space.The algorithm incorporates a variable number of sub-seafloor resistivity layers by adding and deleting interfaces (allowing between 1 and 11 layers) required to minimize the data misfit.It utilizes the Metropolis-Hastings-Green acceptance criterion (Metropolis et al., 1953) to accept or reject the new model state.Additionally, model parameters are perturbed around the current parameter by proposing from a Cauchy distribution, and accepted according to the Metropolis-Hastings criterion.To improve the starting model for MCMC sampling, we use a non-linear optimization technique (Dosso et al., 2001).
The PPD relies on the selection of appropriate data errors and minimum errors.If the data errors are large, the model resolution will decline, leading to large uncertainties.Conversely, small data errors can cause overfitting, resulting in unrealistically low uncertainties.The uncertainty arising from data stacking of the processed data is smaller than systematic measurement errors, which can be due to deviations in array geometry, timing discrepancies between transmitter and receiver signals and positioning of the dipoles, etc.To account for these systematic errors, we apply a minimum relative error of one percent.The one percent value is close to the data error derived during the CSEM data processing.It is important to note that the applied minimum errors are smaller than the minimum allowed relative errors of four percent that were assumed for the deterministic 2-D inversion in Micallef et al. (2020;Figure 2a).This larger error was employed in the 2-D inversion as it searches for a vertically and laterally smoothed model over all data collected along the survey lines, whereas our 1-D approach uses only data measured at one predefined waypoint.To better understand the impact of the applied error model on the PPDs and allow for a more coherent comparison with the 2-D inversion, we also provide PPDs for an assumed minimum error of four percent in the appendix (cf., Figure A1).

Probability Density of Resistivity Models
Inherent to CSEM data inversion is that the product of the resistivity and layer thickness is often better resolved than the parameters themselves (Edwards, 1997).Therefore, one major advantage of implementing a transdimensional inversion is the assessment of the uncertainty of the resistivity and the layer thickness.We present these in terms of interface probabilities as a function of depth (Figure 2b, left) and PPD profiles for a layered sub-seafloor resistivity (Figure 2b, right).To represent the Bayesian results as interface-depth probabilities, we create depth grids and calculate normalized histograms for the number of interfaces estimated at each depth interval (Figure 2b, left).To illustrate the PPDs, we generate normalized histograms for resistivity values at predefined 5 m depth intervals (Figure 2b, right).The histogram values within each depth bin represent the corresponding probability density with higher probability values represented by warmer colors.Cooler and gray colors represent low probability density while white indicates no models within the grid cell.The posterior median model is presented by a dashed black line, and the uncertainties in the model parameters are quantified with 95% credibility intervals (CI).Additionally, to facilitate comparison, we extract resistivity-depth profiles from the 2-D inversion conducted by Micallef et al. (2020) at the respective locations (represented by black circular markers).

Conversion of Resistivity to Pore-Water Salinity Probability Distribution
After obtaining a PPD of the seafloor resistivity through Bayesian inversion, we convert these bulk resistivities into pore-fluid salinity values applying Archie's empirical relationship (Archie, 1942), which links the bulk resistivity (ρ b ) to the pore-fluid resistivity (ρ f ): Here, ϕ is the porosity of sediments, a is the tortuosity factor, and m is the formation cementation factor, indicating the degree of interconnection of the pore spaces and generally assumed to be constant for a given geological structure.In this study, the interpretation of seismic sections (Micallef et al., 2020)  are based on CTD measurements.Within the sediments, we used a local linear gradient of 0.04°C/m, as suggested by in situ measurements at Site U1352 (Fulthorpe et al., 2011).To calculate pressure as a function of latitude and water depth, we use the freely accessible Seawater Library (Morgan, 1994).Subsequently, we compute porewater salinity at different depths based on the estimated pore-water resistivity, temperature and pressure data as implemented by Ruiz-Martinez (2023).
The porosity values in our analysis are derived from two sources: IODP borehole data and seismic facies analysis.At the borehole U1353, porosity values obtained from in situ measurements are available (Fulthorpe et al., 2011).
For locations away from the borehole, we extract the porosity-depth profile at the desired waypoint from the seismically interpreted facies along the survey lines presented by Micallef et al. (2020; Figure 3).We assume that porosity values at each depth are normally distributed, with a standard deviation of 0.05.Our incomplete knowledge of Archie's parameter coefficients a and m is encapsulated in the workflow by assuming a uniform distribution in the range of 0.8-1.2 for a and 1.9 to 2.3 for m.The uniform distribution best captures the variation in Archie's parameters due to interbedding of silty and sandy material, which has been observed in the borehole core and seismic data (cf., Figure 3).We incorporate uncertainties for Archie's parameters, m and a, in our analysis.Both parameters jointly contribute to the overall slope in the relationship between bulk resistivity and pore-fluid resistivity described in Equation 4. The parameter a is theoretically 1.0 as the formation factor F = ρ b / ρ f should be 1.0 at 100% porosity, but has been introduced by Winsauer et al. (1952) to allow a better fit to the sample data and accounts for additional uncertainties (Schwalenberg et al., 2020).Glover (2016) points to measurements errors in porosity, pore fluid salinity, and temperature, which may cause a non-unity value of the parameter a to fit the electrical data well.
To compute pore-water salinity distributions as a function of depth, we generate a depth grid with a vertical resolution of 5 m.For each depth bin, we implement the following steps: 1. Randomly select a bulk resistivity (ρ b ) value according to the resistivity probability function.
3. Randomly select m and a values from their individual distributions.4. Apply Equation 4 to calculate ρ f . 5. Convert the resulting ρ f at each grid to pore-water salinity, following the procedure described previously.6. Compute normalized histograms of the salinity distributions at each depth interval and illustrate them as PDF.
Our classification scheme characterizes pore-water as "fresh" for salinities below 1 psu or equivalent pore-water resistivity above 5 Ωm, "fresh to brackish" for salinities ranging from 1 to 10 psu or pore-water resistivity between 0.5 and 5 Ωm, and "brackish" for salinities ranging from 10 to 30 psu or pore-water resistivity within 0.3-0.5 Ωm.

Workflow Calibration
We apply the previously-described workflow to three key waypoints of the survey area.At each waypoint, the inversion is run until it accepted a range of 80,000-100,000 models which fit the data.We first execute a Bayesian inversion at the waypoint located nearest to borehole U1353 (WP9 on line 2) to verify, validate and calibrate our approach, where we have in situ control.We derive the PPD of resistivity at WP9 as illustrated in Figure 2b, and transform it to pore-water salinity estimates (cf., Figure 4) using (a) the porosity distribution over depth obtained from the porosity-depth profile collected at the borehole and subsequently compare it to (b) the porosity estimates derived from the interpreted seismic facies (cf., Figure 3c).
Subsequently, we repeat the procedure at WP12 and WP60 along line 4, where the most prominent resistivity anomaly is detected, but no borehole data are available to constrain the pore-water salinity conversion (cf., Figures 5 and 6).

CSEM-Derived Salinity Conversion at Borehole U1353
Figure 2a shows the 2-D inversion results of line 2 presented by Micallef et al. (2020).The interface probability (left panel in Figure 2b) at WP9 (line 2) indicates that the CSEM data can resolve two layers over a half-space with interfaces most likely occurring at 14 and 80 mbsf.In Figure 2b (right panel), the upper boundary of a subtle resistivity increase is well resolved at approximately 14 mbsf, while the lower boundary (at approximately 80 mbsf) is less decisive, which is typical for CSEM resolution characteristics.The resistivity PPD (right panel in Figure 2b) denotes that the resistivity of the upper two layers is well-constrained at 0.8 ± 0.5 and 2 ± 0.8 Ωm, respectively, down to approximately 80 mbsf.However, at greater depth, the CI are wider, implying that the CSEM data do not constrain the model resistivity beneath the second layer.The resistivity profile derived by the smoothed deterministic 2-D inversion at this location also shows changes in vertical resistivity-depth gradient.
For the second layer, the 2-D resistivity model at WP9 coincides with the resistivity of the highest probability in the Bayesian inversion.For the layer above and the region below the second layer, the 2-D model predicts higher resistivities than the resistivities at the center of the probability distribution, yet maintains the 95% credibility interval of the PPD.
The pore-water salinity probability density computed for porosity measurements taken in borehole U1353 (Figure 4a) shows a salinity-depth variation with decreased salinity values within the second layer.Above layer two, the predicted salinities approach normal seawater values.Below layer two, the range of possible salinities widens and ranges from brackish to hypersaline (>30 psu) water.Yet, the highest probabilities are shown for seawater salinity of around 34 psu.The median pore-water salinity model follows the salinity trend over depth measured at the borehole (pink circular markers), which indicates a salinity decrease from 34 to 24 psu between 20 and 60 mbsf.The salinity profile converted from the 2-D inversion model fits well to the center of the salinity PPD for both borehole-based (Figure 4a) and seismically derived porosity estimates (Figure 4b) within the aquiferous layer.Salinity predictions based on porosity distribution extracted from interpreted seismic facies (Figure 4b) show almost identical results, validating that seismically derived porosities distributions are applicable to capture the in situ porosity profile, particularly if we assume a Gaussian porosity distribution at each depth for seismic-derived porosities.

Characteristics of Salinity Anomalies at the Center of the Survey Area
The PPD of resistivity at WP12 in Figure 5b shows two well-identified interfaces at depths of approximately 20 and 40 mbsf, each associated with a distinct increase in resistivity.The resistivity within the third layer (ranging from 40 to 150 mbsf) reaches the maximum prior resistivity of around 200 Ωm.Although resistivities below 200   D inversion model in Figure 5a shows only a continuous smooth increase in resistivity with depth.While resistivities of the upper two layers and lower half-space coincide with the median resistivity of the 1-D inversion results, the 2-D inversion model does not include the highly resistive intermediate layer.We attribute the difference in resolution to the lower minimum relative error of one percent assumed in the Bayesian inversion compared to the minimum relative error of four percent applied in the deterministic 2-D inversion.This hypothesis is corroborated by the fact that a rerun of the Bayesian inversion at WP12 with an error level of four percent leads to a simple two-layer model consisting of a conductive layer over a resistive half-space with no resistive anomaly (cf., Figure A1a).
At WP60, the PPD of resistivity (Figure 5c) exhibits three layers over a half-space.The Bayesian inversion places a resistive layer between 30 and 75 mbsf (with the highest probability for a resistivity value of 15 Ωm) over a conductive layer (with resistivity values of less than 1 Ωm) between 75 and 100 mbsf.Note that the resistivity value within the underlying half-space is poorly constrained, allowing for values greater than 3 Ωm.The 2-D inversion model at this waypoint shows a smoothed version of the Bayesian-derived model, including the more conductive layer at around 75 m depth.
Overall, we observe a decrease in resistivity values within our Bayesian resistivity distributions from the near coastal site WP12 to the farther coastal site WP60.
The pore-water salinity distribution derived from the resistivity PPDs at WP12 (Figure 6a) indicates that the salinity values near the seafloor are around 34 psu, which is typical of the seawater salinity of the survey area.A layer of very low salinity appears between approximately 40 and 150 mbsf, with a thickness of 110 m and salinity variations between 0.1 and 0.5 psu indicating freshened pore-water.This layer of fresh pore-water roughly corresponds to the vertical extent of low-porosity (35%) coarse-grained sand deposits in the interpreted seismic section (orange zone in Figure 3b).The correlation of the resistive anomaly with coarse sand structure further raises confidence in our results and that the CSEM data can resolve a layer with freshened pore-water when using an error level of one percent.At a depth between 110 and 200 mbsf, salinities increase (between 0.1 and 10 psu) indicating fresh to brackish pore-water, which extends to greater depths.
At WP60, there is a seawater-saturated shallow layer that is approximately 20 m thick, which is underlain by a low salinity layer extending down to about 75 mbsf.Pore-water salinity in this layer is concentrated around 3 psu, suggesting the presence of a freshened zone between 25 and 75 mbsf (Figure 6b).Note that the depth interval of the freshwater zone occurs at a similar depth interval of the freshened groundwater body identified in the borehole on line 2. The freshened layer is followed by a highly saline layer between 75 and 100 mbsf, with salinity values ranging between approximately 30 and 100 psu.Pore-water salinity decreases again at depths greater than 100 mbsf, showing a variation between 0.3 and 30 psu, suggesting a fresh to brackish groundwater body.These values are slightly higher than the predicted salinities within the lower freshwater body identified beneath WP12 at depths greater than 110 mbsf (Figure 6a).

Discussion
Constraining salinity distributions of pore-fluid using CSEM data requires an understanding of the uncertainty in electrical resistivity, as well as a assuming distributions for the host-rock porosity and Archie's parameters.Deterministic inversion approaches fail to provide a quantitative uncertainty estimate of the resistivity.In contrast, Bayesian inversion algorithms result in the PPD of resistivity, from which uncertainty in pore-water salinity estimates can be addressed.In this study, we present a workflow that converts sub-seafloor resistivity models to a PPD of pore-water salinity by incorporating seismic attribute data and considering uncertainties in sediment porosity and Archie's parameters in Equation 4.

Comparing 2-D Deterministic and 1-D Bayesian Inversion Results
At the borehole location (WP9 of line 2), the probability distribution of resistivity and 2-D resistivity model coincide very well (Figure 2b), which validates our 1-D Bayesian approach.Similarly, at the eastern WP60 on line 4 (Figure 5c), where the deterministic 2-D model exhibits lateral resistivity variations, the PPD of resistivity shows good congruence with the 2-D inversion model at WP60 indicating that the 1-D Bayesian approach can work in regions of moderate 2-D model variations.Surprisingly, at the western part of line 4 (WP12), where the sub-seafloor appears to be 1-D, there are clear disparities between the deterministic 2-D resistivity model and the
probability distribution of resistivity between approximately 40 and 150 mbsf (Figure 5b).The differences can be reconciled when applying a consistent minimum relative error of four percent in the Bayesian and 2-D inversion (cf., Figure A1a).The resistivity model obtained from the 2-D inversion and the PPD of resistivity at WP12 in Figure A1 show similar behavior, although the 2-D resistivity model indicates lower values as compared to the median resistivity model of obtained through Bayesian inversion.This shows that deterministic inversion approaches do not resolve sudden changes in resistivity due to the smoothing regularization of the method.
The comparison between the probability distribution of resistivity depicted in Figure 5 with a minimum relative error of one percent and Figure A1 with a minimum relative error of four percent reveals that the PPD of resistivity is highly dependent on the choice of error.Furthermore, increasing the data error leads to a decrease in model resolution.While we have verified the implementation of the one percent minimum relative error on data at the borehole location, it cannot be ascertained that we are not overfitting systematic errors in our data by choosing a small error of one percent compared to error estimates from stacking during CSEM data processing.The local error at one waypoint is usually smaller than the overall assumed error on the entire study area.Therefore, assuming a smaller minimum relative error and applying the Bayesian approach improves the model resolution.
The difference between the probability distribution of resistivity and the 2-D model underlines the need for conducting additional inversions at more waypoints identified by the 2-D inversion results.Nevertheless, a significant limitation of the Bayesian inversion approach is the high computational time required to calculate the PPD.Consequently, this approach is currently only suitable for a few investigation points and is dependent on a multi-dimensional deterministic inversion to guide the choice of points for stochastic inversions.

Statistical Derivation of Pore-Water Salinity From PPD of Resistivity
Estimating the pore-water salinity from the probability distribution of sub-seafloor resistivity is challenging due to uncertainties in seafloor resistivity, porosity and Archie's coefficients.Borehole U1353 comprises lithological units that include interlayered sand horizons and silty sand material, each characterized by different Archie's parameters (a and m in Equation 4).To address this issue, we choose a uniform probability distribution for each of these parameters to convert the PPD of resistivity to the probability distribution of pore-water salinity.This approach results in broader salinity probability distributions (cf., Figure 4), allowing us to identify potential uncertainties due to our lack of detailed in situ knowledge to calibrate Archie's parameters.
The salinity-depth measurements at borehole U1353 (indicated by pink circular markers in Figure 4) align with the high probability range of the predicted pore-water salinity distribution, although they do not coincide with the median salinity model.Specifically, salinities between 20 and 60 mbsf are higher compared to the salinity probability distribution's center, while salinities in the region below 120 mbsf are lower due to the presence of a different lithological unit (cf., Expedition 317, Scientists, 2011).These results demonstrate that the uncertainty in Archie's parameters and porosity impact the pore-water salinity estimation.

Correlation Between Pore-Water Salinity Distribution and Seismic Facies
The agreement between the probability density distributions of pore-water salinity obtained from resistivity models employing seismically-derived porosity estimates and the in situ porosity measurements at borehole U1353 (Figures 4a and 4b) supports the suitability of using interpreted seismic facies for porosity estimation at distant locations from the borehole.It is important to note that the seismic facies classification needs to be extrapolated from a borehole and may lose accuracy at larger distances to the borehole.
The pore-water salinity PPD at the borehole location (WP9) overlaid on MCS line 2 demonstrates that the first and second interfaces at 14 and 80 mbsf correspond with two seismic reflectors at the same depths (cf., Figure 7a).This alignment suggests that the salinity decrease occurs within an onlap sedimentary package that comprises fine-grained sand, which is both over-and underlain by fine silty sand material.Toward the shore, the onlap sedimentary package pinches out, making it unlikely that there is a connection or feeding mechanism from current onshore groundwater bodies at U1353.In contrast, the sandy deposit extends horizontally to the east.However, at borehole U1354, which is located approximately 8 km seaward of U1353 along line 2, the subtle salinity variations are not detectable (cf., Figure 2c in Micallef et al. (2020)).Overall, the geometry of the sedimentary unit and estimated salinity variations suggest that the freshened OFG body along line 2 is most likely a patchy remnant of a paleo groundwater body that formed during previous glacial low stands (U1353 is located landward of the paleo-shelf during the last glacial lowstand).This remnant body may currently experience salinization due to saltwater diffusion from the seaward side.
The 2-D resistivity model for line 4 presents salinity anomalies with greater variations and deeper depths compared to line 2 (cf., Figures 2a and 5a).Pore-water salinity PPD profile obtained from the eastern part of line 4 (WP60) reveals the occurrence of a shallow freshwater zone between 15 and 70 mbsf (cf., Figure 6b).This OFG body is located at a depth similar to the one observed at U1353, but with significantly lower salinity values.The OFG body is situated within sandy sediments, which are underlain and overlain by silty and fine-grained material.
As with line 2, the seismic reflector that corresponds to the upper limit of the sediment package pinches out landwards.This observation is further supported by the 2-D resistivity section, which shows the extension of the resistivity anomaly from approximately 24 km from the beginning of line 4 to the end of the profile (Figure 5a).We interpret this feature as a remnant OFG body, consistent with the findings and hydrogeological modeling suggested in Micallef et al. (2020), which likely formed during the last (and possibly preceding) glacial seawater lowstand, similar to line 2. A possible geological explanation for the zone of high conductivity (between 75 and 100 mbsf) is the presence of a clay-rich sediment layer consisting of very fine-grained particles located directly beneath the shallow OFG.However, we also acknowledge the possibility of inversion artifacts resulting from 2-D effects, which could provide an alternative explanation.
At greater depths (>100 mbsf) beneath WP60, a second low salinity layer has been identified within silt/fine sand facies (Figure 6b).The upper limit of this OFG body coincides with a seismic reflector located at approximately 100 mbsf, which is present throughout the entire profile and rises landward (cf., Figure 7b).The Bayesian inversion results suggest that the shallow and deep OFG bodies are not interconnected.This is further supported by both the 2-D inversion results and the seismic data, corroborating the absence of hydraulic continuity between these two distinct OFG bodies.The resolution of CSEM data decreases with depth and is limited below 100 mbsf, and as a result, we are unable to accurately constrain the depth of the base of the OFG body and variations in salinity at greater depths.Yet, our results indicate that pore-water salinities in this region are below 10 psu.Water FAGHIH ET AL.
The deeper low-salinity body beneath line 4 represents a somewhat different scenario from the shallower OFG bodies that are geometrically constrained and associated with onlap sediment packages.At the western part of line 4 (WP12), a freshwater OFG body with salinities below 0.5 psu is predicted between approximately 40 to 150 mbsf (cf., Figure 6a).This zone extends slightly above the top of a coarse sand sediment package (orange zone in Figure 3b), and is associated with a prominent seismic reflector that slopes eastwards along the profile, indicating the top of the deeper OFG body beneath WP60 (cf., Figure 7b).The base of the freshwater body is aligned with the base of the coarse-grained sand sediments at a depth of approximately 150 mbsf.Below the coarser sand zone, the pore-water salinity slightly increases while remaining in the freshwater range below 10 psu.

Interpreting an Extensive OFG Body
The correlation between the observed pore-water salinity variations at WP12 and WP60 ranging from 0.1 to 10 psu strongly suggests a hydraulic connection between the lower OFG body beneath WP60 and the freshwater OFG body to the west at WP12 (cf., Figure 7b).Our Bayesian inversion models, 2-D inversion results, and seismic data all support the possibility of a continuous freshwater body that may be hydrologically linked to a land-based aquifer.While CSEM data is missing across the shallow water coastal transition zone, which would have provided continuous coverage to confirm the land-sea groundwater connection, an EM/GPR study conducted along the Ashburton coast (north of line 4) by Weymer et al. (2020) suggests that the onshore-offshore groundwater system is likely connected.However, the hydrogeological models from Micallef et al. (2020) indicate that the portion of recharge from onshore aquifers is small, and the majority of OFG was emplaced during the last glacial periods.The onshore aquifer is situated within gravel layers of sediment sequences (Dommisse, 2006) that were transported by high-energy braided rivers during the last glacial period (Rowan et al., 2012).Unconnected sand and silt layers within these sequences serve as aquitards (Browne & Naish, 2003).
Given that line 4 is located in an area with modern estuaries of various rivers (cf., Micallef et al., 2020), which are likely remnants of the formerly high-energy braided rivers, and considering the seaward regional groundwater flow, it is plausible that the deep OFG body observed at line 4 is a continuation of the onshore aquifer extending seaward.Furthermore, onshore borehole information in the survey area (cf., Davey, 2004Davey, , 2006) ) provides approximate depths of the onshore multi-layer aquifers between 0 and 50 m, 50-90 m, and greater than 90 m.These depth ranges show a reasonable agreement with the depths of the OFG body along line 4 which strengthens the possibility of an onshore-offshore groundwater connection.By estimating the uncertainties in depth and porewater salinity of OFGs, we can offer additional data inputs and constraints for model parameters, particularly depth and salinity, in hydrogeological modeling.This enhancing the reliability and adaptability of models, enabling them to address a broader range of hydrogeological scenarios.

The Impact of Porosity Variation on the Pore-Water Salinity Estimates
Porosity estimates for a sequence of silt and fine sand layers along line 4 based on interpreted seismic facies are higher compared to those obtained from line 2 as there is a lack of borehole data for ground-truthing.To account for the uncertainty in porosity values, we integrate them into the pore-water salinity estimations by assuming normal distributions for seismically-driven porosity with depth.Yet, porosity estimates may vary beyond this range, especially at lithological transitions (e.g., at WP12 between approximately 50 and 150 mbsf in Figures 3b  and 3e).A lower porosity associated with coarse sand units could result in lower salinities.To validate that the very low salinities estimated are not due to low assumed porosity, we repeated the conversion of seafloor resistivity to pore-water salinity by reducing the porosity for the coarse sand facies by 15% (cf., Figure 8).The resulting PPD of pore-water salinity indicates slightly higher salinity values within the depth interval of the coarse sand layer (0.1-3 psu), yet still falls within the freshwater region.This is in agreement with Micallef et al. (2020), who concluded that the resistivity anomalies derived from 2-D inversion cannot be explained by porosity variations alone, but are indicative of freshened pore-water in the respective sediment sections.

Conclusions
We present a methodology that employs a Bayesian workflow to evaluate uncertainty in pore-water salinity predictions within the Canterbury Bight.Our approach utilizes resistivity estimates obtained from CSEM data, in situ porosity measurements, seismic-facies derived porosities, and Archie's rock physics relationship to characterize OFGs.We utilize a trans-dimensional MCMC algorithm to estimate a probability distribution of resistivity models and interface depth probabilities.By implementing Archie's relationship, we derive the porewater salinity distribution and associated uncertainties from the resistivity models.
The efficacy of our workflow was successfully validated using in situ salinity and porosity measurements at the IODP borehole location, where it accurately captured a zone displaying a subtle drop in pore-water salinity values recorded in borehole U1353.Furthermore, our study reveals that interpreted seismic facies can serve as reliable proxies for classifying porosity estimates extrapolated from nearby boreholes.
In the center of the survey area, the correlation between the PPD of pore-water salinity and the seismic reflection profile along coast-perpendicular line 4 suggests the presence of an OFG body containing brackish to freshened pore-water stored within silty/fine-grained sediments at depths greater than 100 mbsf.Shoreward, the OFG body extends and transitions into facies comprising coarse sand sediments (west of line 4 at depths below 40 mbsf, Figure 7b).This freshened OFG likely represents an extension of the onshore aquifer and appears to be disconnected from a shallower local freshened zone identified in the eastern part of the profile.
To enhance our understanding of the extent of the salinity anomalies associated with freshened pore-water, we integrate these anomalies with derived from reflection seismic data.The trans-dimensional approach allows for the identification of abrupt resistivity changes related to freshened pore-water layers which helped to identify facies possibly connected to onshore groundwater that a 2-D deterministic model could not provide.Therefore, we strongly recommend assessing uncertainties in resistivities using Bayesian approaches for selected waypoints and estimating uncertainties for derived salinities to avoid both over and under-interpretation of the hydrogeological model such as the extent of groundwater within the seafloor.Our methodology offers a framework for estimating the salinity and thickness of OFG bodies and assessing their associated uncertainties in the salinity predictions to accurately evaluate volumes of freshened pore-water beneath the seafloor.

Appendix A: PPD of Resistivity for a Minimum Relative Error of Four Percent
To assess the impact of error models on the PPD of resistivity obtained through Bayesian inversion, we perform additional Bayesian inversions at WP12 and WP60 using a minimum relative error of four percent, as applied by Micallef et al. (2020).Figure A1 illustrates the Bayesian inversion results at WP12 and WP60 in terms of interface-depth probabilities and the marginal probability density profile of the resistivity.The probability distribution of resistivity at WP12 only resolves the first interface at 25 mbsf, indicating a significant increase in resistivity to values greater than 100 Ωm.Overall, the data lose resolution at greater depths.At WP60, the PPD of resistivity shows three distinct layers above the basement, which agrees well with the 2-D inversion model.Nevertheless, interface probabilities as a function of depth suggest only one distinct interface at 50 mbsf with high uncertainty.

Figure 1 .
Figure 1.(a) Map of the Canterbury Basin displaying the locations of acquired CSEM and seismic profiles.The yellow circles illustrate the locations of the IODP boreholes.CSEM line 2 intersects with the location of IODP site U1353, for which (b) salinity and (c) porosity measurements up to a depth of 300 mbsf are available (Fulthorpe et al., 2011).

Figure 2 .
Figure 2. (a) 2-D bulk resistivity model derived from acquired CSEM data along line 2 as overlain on the corresponding seismic reflection profile (Micallef et al., 2020).Depth is shown in meters below sea-level (mbsl) with a vertical exaggeration of 22. Black triangles indicate stationary waypoints for CSEM data collection.The location of borehole U1353 and the closest waypoint (WP9) to that is marked by an arrow.(b) The result of 1-D Bayesian inversion at WP9.The left panel shows interface probability as a function of depth.The right panel presents the resistivity marginal probability profile.The color indicates the probability.Credibility intervals contain 95% of the model samples evaluated at each depth interval.The extracted 2-D bulk resistivity model is shown by black circular markers.Depth in this panel is in meter below seafloor (mbsf) and horizontal gray dashed lines indicate the depth of interfaces.(c) Stacked inline electric field data measured at WP9 with gray error bars at four receivers (Rx1-Rx4).The forward responses of 20 randomly selected models from the resistivity model distribution are presented in blue.

Figure 3 .
Figure 3. Seismically interpreted facies along (a) line 2 and (b) line 4 from Micallef et al. (2020) with a vertical exaggeration of 22. (c) Porosity profile measured at borehole U1353 (blue circular markers) and porosity profile extracted from interpreted facies (yellow circular markers) at the borehole location.The assumed probability distribution of porosities over depth at (d) WP9, (e) WP12, and (f) WP60.

Figure 4 .Figure 5 .
Figure 4. (a) Distribution of pore-water salinity as a function of depth at WP9 derived from Bayesian bulk resistivity probability density (Figure 2b) by applying Archie's relationship using (a) porosities derived from in situ measurements at U1353 and (b) porosities derived from interpreted seismic facies.Pore-water salinity measurements at borehole U1353 and salinity estimates converted from 2-D inversion bulk resistivity models are shown in pink and black circular markers, respectively.Horizontal gray dashed lines represent the depths of resolved interfaces.

Figure 6 .
Figure 6.Probability distribution of pore-water salinity converted from posterior probability density (PPD) at (a) WP12 and (b) WP60.The pore-water salinity distributions are derived based on the PPD of bulk resistivity (Figures5b and 5c) and interpreted seismic facies derived porosity distributions at the respective waypoints (Figures3e and 3f).Horizontal gray dashed lines in both panels represent the depths of resolved interfaces.
Figure 6.Probability distribution of pore-water salinity converted from posterior probability density (PPD) at (a) WP12 and (b) WP60.The pore-water salinity distributions are derived based on the PPD of bulk resistivity (Figures5b and 5c) and interpreted seismic facies derived porosity distributions at the respective waypoints (Figures3e and 3f).Horizontal gray dashed lines in both panels represent the depths of resolved interfaces.

Figure 7 .
Figure 7. Multi-channel seismic reflection data along (a) line 2 and (b) line 4 showing seismic reflectors.The white line in (b) shows the boundary between two seismic profiles along line 4 acquired using different acquisition equipment and geometry.Probability density distribution of pore-water salinity at (a) WP9 and (b) WP12 and WP60 overlain on the corresponding seismic lines.

Figure 8 .
Figure 8.(a) An assumed 15% decrease in porosity within the depth range of coarse sand sediments at WP12 along line 4 (indicated by the orange zone in Figure3b) between approximately 50 and 150 mbsf.The corresponding probability density distribution of pore-water salinity is calculated for a porosity value of (b) 35% and (c) 20% in the coarse-grained materials.Horizontal gray dashed lines present the depths of interfaces.

Figure A1 .
Figure A1.Probability density distributions of sub-seafloor resistivity at (a) WP12 and (b) WP60 assuming a minimum relative error of four percent, as it was assumed for 2-D deterministic inversion of CSEM data at the same waypoints.The depths of interfaces at each waypoint are highlighted by horizontal gray dashed lines.