An integrated modeling scheme for characterizing 3D hydrogeological heterogeneity of the New Jersey shelf

Continental shelves around the globe are hosts to vast reservoirs of offshore freshened groundwater. These systems show considerable complexity, often as a function of the geological heterogeneity. Data needed to characterise these systems are often sparse, and numerical models rely on generalized simplifications of the geological environment. In order to improve our understanding of these systems, it is necessary to implement modelling approaches that can produce large-scale geologically representative models using sparse data. We present an interdisciplinary stochastic modelling workflow incorporating borehole data, 2D depth-migrated seismic profiles, seismic attributes, and prior knowledge of the depositional setting. We generate a conditioned Gaussian field of porosity on the New Jersey shelf. We also perform a petrophysical conversion to a corresponding permeability distribution. The model dimensions are 134 km × 69 km × 1.7 km, with an adjustable resolution that can be adapted for process-based models of flow and solute transport. The integrated approach successfully translates small-scale porosity variations to a shelf-scale model that honors key characteristics of the New Jersey shelf wave-dominated depositional environment. The model was generated using open-source packages. All data and code to reproduce the complete workflow are provided along with this study so the model can be reproduced at any resolution for further studies of continental shelf processes offshore New Jersey.


Introduction
Offshore freshened groundwater (OFG) refers to fluids stored in sediments below the seafloor that are less saline than ocean water. This phenomenon has been observed in several continental margins around the globe (Post et al. 2013). The documented occurrences show that these systems vary considerably in their depth range, geometry and distance from the shoreline (Micallef et al. 2020a). The characterization of offshore freshened groundwater on continental shelves requires integration of multiple subsurface data sets. Several studies have shown that hydraulic properties of the subsurface are strongly influenced by the geological heterogeneity (Simmons et al. 2001;Michael et al. 2016;van Engelen et al. 2019;Zamrsky et al. 2020). Numerical forward model simulations have become one of the key tools in understanding the emplacement and subsequent survival of offshore freshened groundwater systems. As numerical modeling techniques have improved, increasingly complex flow models have been used to investigate OFG systems in various environments (Cohen et al. 2010;Siegel et al. 2014;Thomas et al. 2019;Micallef et al. 2020). A key factor in any reliable model is capturing the spatial variation of flow properties on a meaningful scale.
The correlation between subsurface geological features and freshened pore-water has been observed across various depositional environments. Grasby et al. (2009) analyzed pore fluids in the Beaufort-Mackenzie basin in Ariel T. Thomas and Jan von Harten contributed equally to this work.
11 Page 2 of 19 northern Canada, a post-rift infill of sediments deposited in a marine environment. It was observed that regional faults delineate the spatial distribution patterns of pore fluid salinity. Numerical simulations on continental shelf models have demonstrated the impact of low-permeability sediments on the salinification of offshore reservoirs in Suriname (Groen et al. 2000) and New Jersey (Cohen et al. 2010;Siegel et al. 2014;Thomas et al. 2019). A recent study by Micallef et al. (2020) discovered that significant lateral variability of formation resistivity was correlated to shore-parallel heterogeneity in subsurface properties. This highlights the importance of constraining heterogeneity in three dimensions for numerical forward modeling. However, data availability as well as expensive computational requirements for complex 3D numerical simulation are often limiting factors. Therefore, 2D models have been widely used in the investigation of offshore fresh groundwater systems, some examples include the New Jersey shelf (Thomas et al. 2019), offshore Massachusetts (Person et al. 2003;Marksamer et al. 2007) and Greenland's continental shelf (DeFoor et al. 2011). These studies have shown that 2D models can provide meaningful insights on the dynamic nature of fresh and saline water interaction. However, 3D models can account for additional factors influencing freshwater distribution, such as seafloor topography, as highlighted by Siegel et al. (2014) in a numerical study of OFG emplacement offshore Massachusetts during glacial periods. The increase in computational capability and efficiency of numerical algorithms for variable-density flow, such as distributed memory parallel computing (Verkaik et al. 2021), allows for the consideration of realistic 3D heterogeneity in subsurface models.
Geological modeling often faces the challenge of extrapolating subsurface properties from sparsely distributed data. Seismic reflection data has been widely used to overcome this challenge and to achieve reliable geological models. Seismic reflection data play important roles, not only for defining the structure and stratigraphy of the models (Shtivelman and Goldman 2000; Pugin et al. 2009;Rasmussen 2009;Micallef et al. 2020), but also to estimate the fluid contents as well as the distribution of rock-physical properties such as density and porosity (Batzle and Wang 1992;Al-Bazali et al. 2008;Bosch et al. 2010;Jusri et al. 2018). A common way to utilize post-stack seismic reflection data to determine rock-physical properties is through a seismic acoustic impedance inversion (Latimer et al. 2000;Yilmaz 2001). Furthermore, the acoustic impedance can be integrated in a seismic multi-attribute analysis to predict subsurface rock-physical properties, such as porosity, with a higher degree of accuracy (Hampson et al. 2001). A previous study has shown that the same method is applicable in estimating porosity distribution of fresh water reservoir at the New Jersey Shelf with promising results as well (Jusri et al. 2018).
This study presents an innovative workflow for characterizing the 3D heterogeneity of flow properties on the New Jersey shelf, which integrates the key characteristics of the wave-dominated depositional setting. Seismic interpretation of key sequence stratigraphic units on 2D depthmigrated profiles defines the structural framework of the model. The geo-structural model was generated using the open-source geological modeling toolkit GemPy (de La Varga et al. 2019). We integrated borehole data and seismic attributes in a geostatistical modeling workflow to determine the porosity in the model domain outlined in Fig. 1. We incorporated the seismic attribute as secondary porosity data using a pixel-based analysis to generate a spatially extensive data-set for variogram modeling. The porosity distribution was generated as a Gaussian random field (Heße et al. 2014), using the open-source toolkit GSTools (Müller and Schüler 2021). Petrophysical relationships derived from laboratory measurements of porosity and permeability presented in Lofi et al. (2013) were then used to translate the simulated property model to a permeability distribution. This study presents a parameterized model of the New Jersey shelf that captures realistic 3D geological heterogeneity that can be used to study continental shelf processes such as paleo-hydrogeological reconstruction of fresh groundwater on the shelf, or modeling of the shelf geomorphology. The workflow and resources provided within this study can be used to easily reproduce stochastic model realizations for future studies.

Depositional setting and heterogeneity
The New Jersey margin has been characterized as a wavedominated depositional environment (Pekar et al. 2003). The New Jersey shelf clinothem complex was more recently described by (Proust et al. 2018) as a compounded depositional system in which a sub-aerial delta feeds a more distal subaqeous counterpart, via a sediment bypass controlled primarily by the fair-weather wave base. In a wave-dominated setting, wave processes drive the redistribution of most of the sediments supplied to the delta front. Therefore, progradation involves the entire delta front, not only concentrated points. This forms laterally amalgamated facies, which are near linear and sheet-like features parallel to the shoreline (Reading 1978;Willis et al. 2021). (Monteverde et al. 2008 presented isopach maps of 20 Miocene sequences on the New Jersey margin based on 2D marine seismic profiles. The isopachs show evidence of nearly linear shoreline-parallel depositional centres indicative of along-shore sediment reworking by wave activity. (Monteverde et al. 2008) describes this as being consistent with the wave-dominated margin model.
Wave-dominated environments are known in the petroleum industry to be associated with high lateral reservoir continuity and internal heterogeneity along the strike direction (Sech et al. 2009;Olsen et al. 2017;Berton et al. 2021). On the New Jersey shelf this is evidenced by the parallel to sub-parallel seismic reflector patterns observed on strike lines, consistent with sheet-like deposition across the model area. Considering this depositional environment, we postulate that the porosity distribution would also display similar geometric trend, with longer correlation lengths in the strike direction (along-shore) relative to the dip direction.

Methodology
This section describes the integrated approach to generate a representative 3D hydrogeological model of the New Jersey shelf. Firstly, we show the seismic interpretation of key sequence boundaries, which forms the structural framework of the model. This is followed by the geostatistical analysis of well log and seismic attribute data. Well log data is used to determine the vertical variogram properties, whereas the horizontal variogram properties are derived from seismic attribute data. Finally, we present the details of a Sequential Gaussian Simulation scheme used to simulate porosity distribution within the deterministic stratigraphic framework. The approach is summarized in Fig. 2.

Input Data
The primary porosity data for the model was obtained from three boreholes M27, M28 and M29, drilled in 2009 during International Ocean Discovery Program (IODP) Expedition 313 (Mountain et al. 2020). The porosity refers to a calculated porosity derived from the difference between the wet and dry weights and the dry volume of sediment samples. The samples were at approximately 3 m intervals. In the shallow section of wells M28 and M29, poor recovery of mostly unconsolidated sediment resulted in some missing measurements.
The raw seismic data were obtained from previous seismic surveys, i.e., R/V Oceanus cruise Oc270 acquired in 1995 (Mountain 2008) and the R/V Cape Hatteras cruise CH0698 in 1998 (Mountain et al. 2018). The 2D seismic data were depth migrated using a reflection seismic workflow focused on tomographic velocity model building and pre-stack depth imaging. The interval velocity model used for depth migration was computed by several iterations of geologically-constrained reflection tomography as described in a previous study by (Riedel et al. 2018).

Seismic interpretation and geo-structural modeling
We performed a sequence stratigraphic seismic interpretation on 2D depth migrated seismic profiles. The seismic profiles provide images of the shelf sequences in strike and dip orientations. The interpretation consists of seven key sequence boundaries from the Oligocene through Miocene to the ocean floor. This interval consists of siliciclastic sediments that have been characterized in several drilling campaigns and stratigraphic studies (Miller et al. 1996;Miller and Snyder 1997;Mountain et al. 2020). The sediments rest on a former carbonate platform that was buried by deltaic sediments during a late Cretaceous marine transgression (Browning et al. 2008).
Sequence boundaries on the New Jersey shelf have been correlated with major impedance contrasts that can be identified in the seismic records (Miller et al. 2013). Riedel et al. (2018) presented the workflow for depth migration of the seismic profiles, yielding predicted reflector depths for key sequence boundaries on the line Oc270.529. The predicted depths deviated within 5 m of the true core depth reported in IODP Expedition 313 wells M27, M28 and M29. The interpretation presented in this study builds on these depths, and ties the interpreted seismic horizons across all 2D seismic profiles in the model domain. The structural stratigraphic model is defined by seven key surfaces. The sequence boundaries interpreted across the model domain correspond to the Miocene sequences: m1, m4.1, m5, m5.4, m5.8, m6. The basal surface of the model is an erosional surface, o1, which has been interpreted as the sequence boundary straddling the Eocene-Oligocene boundary (Miller et al. 2013). The interpretations of the dip oriented profiles are shown in The surfaces were generated in Petrel using the convergent interpolation method on the 2D interpreted framework. The surface points were then exported to the opensource Python library GemPy (de La Varga et al. 2019), which uses an implicit modeling approach (Lajaunie et al. 1997;Calcagno et al. 2008) to create a full 3D structural sequence model of the New Jersey shelf. The model domain has dimensions 134 km × 69 km × 1.7 km. The modelled sediments range from a thickness of 200 m at the landward boundary of our model to 1700 m in a seaward thickening wedge. The geostructural model presented in this paper has a discretization of 138 × 268 × 85 cells, for a total of 3.14 million cells, leading to a model resolution of 500 m × 500 m × 20 m. The accompanying code can be used to reproduce

Geostatistical modeling
The geostatistical analysis presented in this study is conducted using the two sources of information available: (a) borehole data from IODP expedition 313 and (b) 2D seismic depth profiles. The porosity data, measured on core samples from three boreholes by Mountain et al. (2020), serves as the basis for vertical variogram modeling, and the input for subsequent conditioned geostatistical simulation. As borehole data is insufficient to estimate horizontal variograms, a seismic attribute will be used to estimate the correlation length in these directions. The geostatistical modeling approach incorporates prior knowledge of the New Jersey shelf stratigraphic setting, which is used to provide soft constraints of the vertical and lateral facies patterns, and qualitatively assess the viability of stochastic realizations (Labourdette et al. 2008). The open-source Python library GSTools was used for the geostatistical modeling (Müller and Schüler 2021).
Stratigraphic units were modelled individually, however, global assumptions had to be made for some geostatistical parameters due to sparse data. A single vertical variogram model was created based on borehole data from M28 for all sequences. M28 was selected because of its position on the shelf penetrating topsets, bottomsets and rollover sediments. We assume this well to be representative of the vertical variability in the deltaic depositional system. A single rotational angle was defined for each of the horizontal orientations in each stratigraphic sequence. As a result, we simplify the curvi-linear geometry of the clinothem. Therefore, our analysis may underestimate the property correlation length in some regions of the shelf.

Data preparation
The porosity log data was de-trended and normal-score transformed prior to variogram modeling, followed by the generation of the conditioned random fields. Two trends were detected: (a) A vertical compaction trend modelled using an exponential regression and (b) a trend of increasing average porosity in the offshore direction modelled as a linear trend. Shelf sediments typically undergo compaction due to the weight of overburden material, which results in a reduction of porosity with depth. The porosity data of wells M27, M28 and M29 show a slight reduction of porosity with depth. The data was fitted to the following exponential trend, where is the porosity at depth z (m), 0 is the initial porosity and a is the correlation coefficient. This type of trend, Fig. 4 Seismic interpretation of sequence boundaries shown on strike sections. These profiles, oriented North-East to South-West, show the continuous, parallel to sub-parallel horizons first described by (Athy 1930), is typically used to represent compaction. The In wave-dominated delta systems like the New Jersey shelf, the facies trends are relatively predictable. Progradation of the shoreface generally results in a basinward transition from coarser to finer grain sediments (Stanley and Swift 1976;Niedoroda et al. 1985). Integrating such sedimentary concepts can help to reduce uncertainty of model predictions away from well locations, and ensure that model realizations are compatible with the depositional environment. This geological information can only be accounted for if it is made quantitative (Marion et al. 2000). In our study, we calculated the average porosity at the three borehole locations and ftted a linear trend to the average porosity as a function of distance offshore. The linear trend is defined as where avg and x represent average porosity and offshore distance, respectively. c and m are the intercept and slope of the linear fit.
Normal-score transformation of the de-trended data was conducted using a quantile-quantile transformation, aimed at producing a normally distributed data set. This transformation is a necessary criterion for conditioning data used in the generation of Gaussian random fields in the subsequent steps (Deutsch and Pyrcz 2014).

Vertical variogram analysis of well data
An experimental variogram for the vertical direction was created using a maximum lag distance of 253 m which is one third of the domain, corresponding to the depth range of borehole data from M28, and a bin size equal to the data spacing of 3 m (Deutsch and Pyrcz 2014). For the theoretical variogram, the three basic models of Gaussian, Spherical and Exponential variograms with a nested nugget effect were compared for determining the best fit. Resulting fits of the theoretical models to the experimental variogram were compared using a standard coefficient of determination score, r 2 .

Lateral correlation length derived from seismic data
The calculation of directional variograms from the primary porosity data in this study is restrained by the limited borehole data. We overcame this problem by establishing a correlation between a relative acoustic-impedance attribute and porosity at the well locations. The attribute is defined in Petrel as the sum of regularly sampled amplitude values calculated by integrating the seismic trace (Schlumberger 2010). Based on this correlation, we make the assumption that lateral changes in the amplitude of the attribute reflect changes in the porosity. Therefore, correlation lengths derived from variogram analysis of a secondary property, the seismic attribute, can be used to model the porosity. This approach allows us to determine a vertical to horizontal anisotropy ratio of the correlation length in each sequence. Directional horizontal variograms were calculated along two key orientations of the shelf system, along-shore (NE-SW) and down-dip (NW-SE). Two representative seismic sections, the dip-oriented line OC270.529 and the cross-line in strike direction, CH068.22 (see Fig. 5), were used for this analysis. These lines were chosen due to their central location in the model domain and proximity to IODP wells M27, M28 and M29. The details and plots of the associated variograms can be examined in the variogram modeling section of the accompanying code. In this paper, we focus on the derived correlation lengths assigned to the stochastic model.
We calculated the Spearman's rank correlation coefficient (Hartmann et al. 2018), between porosity and the value of the pseudo acoustic impedance attribute extracted at the approximate well locations. These values were extracted using digital image analysis of the attribute section. The image was first converted from the Red-Green-Blue (RGB) color space to grayscale, to obtain a single integer representation of the color intensity. This conversion was done using a common approach known as the luminosity method (Güneş et al. 2016), where Grayscale = 0.299 × R + 0.587 × G + 0.114 × B . Both data sets were then normalized between [0,1], a standard practice for quantifying correlation between seismic attributes and well log data (Deutsch and Pyrcz 2014). Figure 6 shows the porosity The variogram analysis of the seismic attribute data was performed on the normalized pixel amplitude. We randomly sampled the data of each sequence individually. An example of the randomly sampled data set is shown for the strike and dip directions of sequence m5.8 in Fig. 7. The sample size was set to 25% of the total number of pixels in the image, resulting in approximately 10,000 samples per sequence. The sample data sets for all stratigraphic sequences are included in the appendix, Fig. A3 in Supplementary Material.
For both examined horizontal directions, a rotational dip angle (compare Fig. 7) of highest continuity was determined per sequence by iteratively computing directional variograms to capture the geometry of the shelf setting. Experimental variograms were calculated using a maximum lag distance of half the section length, a bandwidth of 2 m, a tolerance angle of 22.5 circ and a bin size of 20 m. The assign values represent general tolerance parameters for variogram calculation (Chiles and Delfiner 2012). These parameters generated experimental variograms that showed structure and were consistent with the geologic setting. In order to guarantee robust variograms, only those with at least 2500 point pairs per bin were accepted. For the down-dip (NW-SE) direction, we scanned for values of between 0 circ and 2 circ in steps of 0.01 circ , representing typical dip angles for the shelf sequences (Proust et al. 2018). Along-shore (NE-SW), a variable angle range for between − 0.5 circ and 0.5 circ in steps of 0.01 circ was covered to detect possible cross-stratification features. Lateral migration of depo-centers in deltaic sediments can lead to seismic reflectors dipping in opposite directions. In accordance with the vertical variogram, a theoretical model was fitted to the resulting directional experimental variograms, and for each sequence the highest correlation length and the corresponding rotational angles and were recorded for use in the further geostatistical modeling process.

Random field generation
Conditioned Gaussian random fields were created per sequence using the Randomization method (Heße et al. 2014;Kramer et al. 2007;Dentz et al. 2002), as implemented in GSTools (Müller and Schüler 2021). Theoretical variograms and geometric anisotropies estimated in the previous steps were applied to field generation. In addition, the de-trended and normal-score transformed porosity data from wells M27, M28 and M29 was used as conditioning data to constrain random field realizations. Finally, resulting porosity values of the created Gaussian random fields Fig. 6 Well log plots showing the porosity and synthetic log of pixel amplitude of the relative acoustic impedance attribute at IODP well locations M27, M28 and M29. The true values have been normalized on a scale between (0,1) for the analysis. The Spearman's rank correlation coefficient is shown for each well. The values were -0.41, 0.47 and 0.41 for wells M27, M28 and M29, respectively were back-transformed using the quantile-quantile transformation established in the data preparation step. Removed vertical and horizontal trends were reapplied to the resulting field for a final porosity model (Deutsch and Pyrcz 2014). Random fields were estimated on the grid of the constraining geo-structural model with the aforementioned resolution. In this paper, we will present a single realization of the random function showing the fully condition model. In order to assess the predictive capabilities of the model, we generated a realization in which only two of the three wells were assigned as conditioning data. M28 was chosen as the blind well. This blind test model is compared with the fully conditioned model in the Results. Based on the provided methodology and code, model-ensembles for subsequent work can easily be generated.

Porosity to permeability conversion
The final component of the shelf model is the determination of permeability in the model domain. Permeability contrasts in shelf sediments have been shown to strongly influence the dynamics and distribution of offshore freshened groundwater systems (Micallef et al. 2020a). Permeability measurements in the field are often limited, therefore, several empirical models exist that express permeability as a function of various subsurface properties such as: grain size, mineralogy, porosity and pore dimension (Nelson 1994). In this study, we use the following linear regression to determine a relationship between porosity, , and the logarithm of the permeability, where log(k) is the predicted logarithm of permeability; a and b are the regression coefficients. The linear equation would be used to obtain a corresponding permeability distribution for each stochastic Porosity model realization. In the study area, there are ten permeability measurements, conducted on core samples by , this data is summarized in Table 1. Figure 8 shows the cross-plot of log permeability versus porosity values based on data shown in Table 1. The resulting linear fit is defined log(k) = 3.5 + 9.3 , with an r 2 value of 0.71.

Results
In this section we present the key findings of this study, a geo-structural model, variogram model parameters and a parameterization of flow properties on the New Jersey margin centered on IODP Expedition 313 wells.

Geo-structural Model
The 3D stratigraphic framework of the model is shown in Fig. 9, with a vertical exaggeration of 30. The interpreted interval extends from the present-day seafloor to the sequence boundary o1. This is the deepest characterized sequence boundary that can be interpreted on the seismic data. The bounding surfaces are continuous across the model domain, with no major faults interpreted in this passive Fig. 7 Randomly sampled data points from the relative acoustic impedance attribute cross section for sequence m5.8. The angles and are the angle of highest continuity determined in the variogram analysis for the dip and strike sections, respectively. The dashed grey line shows the intersection of the two seismic lines margin setting. The sequences form a seaward thickening wedge, with the maximum dip-direction being toward the South-East. The 3D model shows that the interpreted sequences are also dipping toward the South-West. However, it is noted that this dip angle is more pronounced deeper in the sediment column. The sea floor topography is relatively flat in this orientation. This can be observed in the cross section shown in Fig. 9c.
The seismic characteristics of the sedimentary units are best observed on seismic line Oc270.529 shown in Fig. 4. The shallow units, m1 and m4.1, are characterized by primarily aggradational stacking pattern. The reflectors are gently dipping in the offshore direction with a hummocky characteristic. The mid-Miocene units, m5-m5.8, display higher angle clinothem systems, consistent with high sedimentation rates and delta progradation in this time period (Miller et al. 2013;Proust et al. 2018). Top-and bottom sets of these stratigraphic units show coherent seismic reflections on the seismic section. The basal units, m6 and o1, display a combination of aggradation with some progradational features toward the landward side of the model. We calculated the true stratigraphic thickness (TST) of each interval, this is shown in Fig. 10. These maps are consistent with wavedominated system, in which sediments were laterally redistributed toward the South-West as the delta prograded across the shelf.

Variogram model
Prior to calculating the experimental variogram, the vertical and horizontal trends were subtracted from the borehole data. The parameters of the vertical trend described by Eq. 1, were 0 = 46.03 and −3.75 × 10 −5 . It represents a slight overall reduction of porosity with depth. The fitted trend line is shown with all the data points in Fig. A4 of the Appendix in Supplementary Material. The horizontal trend calculated was avg = 40.4 + 0.06x . This suggests a small increase of porosity in the offshore direction.
The experimental and theoretical vertical variogram based on de-trended and normal-score transformed primary porosity data from well M28 is shown in Fig. 11. The theoretical model with the best fit is the exponential model with a coefficient of determination score, r 2 , of 0.32 (Gaussian: 0.29, Spherical: 0.23). It features a range of 14 m, reasonably corresponding to the thickness of lithologic units described based on core sample analysis by Mountain et al. (2020).
Results for the iterative estimation of highest horizontal correlation lengths per sequence and key direction based on secondary seismic information are presented in Table 2. A robust variogram could not be estimated for sequence m1, due to a low signal-to-noise ratio of the seismic data. In sequence m5, an insufficient number of point pairs were obtained in the thin interval imaged on the seismic strike line. In these cases we assigned the correlation length of the underlying sequence and set to 0 circ , as no additional information is available.
Down-dip (NW-SE) rotational angles, , vary between 0.07 circ and 0.58 circ , which is in general agreement with the regional depositional setting on a gently dipping platform  (Proust et al. 2018;Miller et al. 2013). While younger sequences feature lower angles and are more horizontal, the mid-Miocene sequences feature more steeply dipping trends consistent with the prograding clinothems. Correlation lengths in this direction range from 1365 m to 3447 m. Along-shore (SW-NE) rotational angles , show small positive values between 0.00 circ and 0.47 circ , meaning a vergence towards SW, for all sequences except m41, which features the longest correlation length in a NE dipping angle of − 0.01 circ . Along-shore correlation lengths are consistently higher than the dip direction, with a range between 3076 m and 7426 m. The anisotropy ratios display consistently higher values in the x direction, i.e., along-shore. Figure 12 shows one realization of the porosity distribution model. The range of porosity values captured in the stochastic model is consistent with the values that have been observed in the field measurements. The reduction of porosity with depth can also be seen, with sequence o1 showing generally lower porosity values. The highest porosity values appear in the middle-Miocene sequences, m5, m5.4 and m5.8. The modelled porosity values are compared to the observed porosity at the well locations in Fig. 13. This plot shows that the model captures a significant amount of heterogeneity in the subsurface flow property, and closely matches the general trend of the measured data.The results of the blind test model are shown as the grey curve in Fig. 13. This shows that both the blind and fully conditioned models capture the overall trend of the porosity at the well locations. The blind model only slightly over-estimates the porosity in M29.

Porosity: Permeability relationship
The linear function was used to transform the porosity model realization (Fig. 12) to the corresponding permeability distribution shown in Fig. 14. The modelled values range from the order of 10 −17 m 2 (or 0.1 mD) at the lowest to 10 −12 m 2 (or 1000 mD) at the high end of the permeability scale. This fits the range of expected permeability values for the New Jersey shelf facies ). Due to the use of a linear transform, the features observed in the permeability distribution mirror the porosity model. While Fig. 14 represents a single realization, some key features of the model are consistent with field observations by previous studies (Mountain et al. 2020;Lofi et al. 2013). This includes zones of low permeability that can be observed just below the seafloor, as well as laterally extended low-permeability features near sequence boundaries.

Discussion
In this study, we present an integrated modeling work flow for the parameterization of shelf-scale flow properties on the New Jersey margin. The workflow incorporates 2D depth migrated seismic data, seismic attributes, geostatistical workflows and petrophysical analysis to arrive at a representative stochastic model, that accounts for geological heterogeneity in the subsurface. The flow model was designed for numerical investigation of the extensive offshore freshened groundwater (OFG) reservoir system that has been identified on the mid-Atlantic margin (Mountain et al. 2020). Geologic heterogeneity, and the associated heterogeneity of subsurface flow properties, has been identified as one of the key factors influencing OFG systems (Yu and Michael 2019;Michael et al. 2016;Zamrsky et al. 2020). Micallef et al. (2020a) discussed the integration of heterogeneity derived from seismic data as a key way to improve 3D characterization of shelf stratigraphy for numerical modeling purposes. Our model achieves this by integrating the geo-structural shelf model, primary  The use of relative acoustic impedance attribute as a secondary data for porosity was designed to compensate for the limited well log data. The seismic data provided data points spread over 70 km in the dip direction and approximately 50 km along strike. Attribute data also allows us to define variogram properties in the along-shore direction, which would otherwise not be possible given the orientation of available wells. It is justified by the strong correlation between chronostratigraphic unit boundaries and impedance contrasts reported by Miller et al. (2013). Impedance contrasts are determined by the product of sediment density and sonic velocity. Porosity significantly influences these seismic properties in the subsurface . Wang et al. (2009) showed that changes in porosity induced similar changes in seismic wave velocities. It follows that the variability of acoustic impedance has underlying dependency on porosity. Acoustic impedance attributes and others, derived from seismic data are often used to estimate subsurface porosity distribution (Hampson et al. 2001;Na'imi et al. 2014;Soleimani et al. 2020). In a stochastic reservoir modeling study presented in (Coburn et al. 2006, p. 253), acoustic impedance is used to derive porosity trend maps for a carbonate reservoir. In that study, absolute values of the correlation coefficients between impedance and porosity Fig. 13 Comparison of fully conditioned model (red), blind test model (grey), and measured porosity at well locations for single stochastic realization ranged from 0.41 to 0.43. The porosity realizations were reported to closely match the original data distributions. Some applications rely on very strong correlation between attributes and well log data to directly invert attribute amplitude to obtain porosity. In our case, we use the attribute as a secondary data to derive variogram properties. Therefore, we tolerate the moderate correlation coefficient values observed between core-sampled porosity values and the relative acoustic impedance attribute (See Fig. 6).

Geo-structural model
The geo-structural model is constrained by seven major sequence boundaries. Two criteria were considered in selecting the sequence boundaries: (1) the sequence boundaries yield coherent seismic reflections that could be easily identified on the seismic sections, (2) they subdivide model domain into sedimentary units with comparable seismic stratigraphic characteristics. Seismic sequence stratigraphy has been widely used to characterize sediments in continental margins (Neal and Abreu 2009;Richard et al. 1985;Vail 1987;Embry and Johannessen 1993). Seismic facies analysis has previously been applied in the study of offshore freshened groundwater systems to estimate reservoir geometry and connectivity Inwood et al. 2013;Micallef et al. 2020).
Our isopach maps of the individual stratigraphic units support the wave-dominated depositional environment of the Oligocene -Miocene deposits on the New Jersey shelf (Proust et al. 2018;Pekar et al. 2003). Understanding the geometry of sedimentary units in the model domain gives insight into the distribution of flow properties across the shelf which is important for the characterization of the OFG system. The South-West trending orientation of the sedimentary units is likely to impact the geometry of freshened groundwater intervals and the dynamics of the system. Mountain et al. (2020) showed that the freshwater offshore New Jersey occurs in a complex multi-layered reservoir. Lofi et al. (2013) showed that most of the freshwater in IODP wells correlated with low-permeability intervals. Gustafson et al. (2019) reported a broad resistivity anomaly on the US Atlantic margin, suggesting a laterally extensive freshened groundwater system. In order to reconcile these observations in a numerical model, an integrated approach must be adopted to account for the characteristic heterogeneity of the system. Our model reflects the key depositional trends and provides a good basis for future numerical investigation.

Geostatistical modeling
Variogram analysis faces two significant challenges: limited well log data and the centralized location of wells relative to the model extent. Geological structures often display locally variable characteristics, such as varying azimuth of anisotropy or dip angles (Boisvert et al. 2009).
Our separation of the model domain into seven sedimentary units, captures some of this variability in anisotropy on a regional scale. Each of the units was constrained by its own set of parameters as summarized in Table 2. By this approach, we integrate the variability in the orientation of sedimentary bedding planes that is characteristic of the depositional environment. (Michael et al. 2016) discussed the importance of capturing geological heterogeneity for modeling offshore freshened groundwater. Their study showed that preferential flow through high-permeability conduits creates complex salinity distributions. Studies have shown that groundwater flow is strongly influenced by the orientation of bedding planes (Michalski and Britton 1997;Fan et al. 2007). Small-scale porosity data, sampled at 3 m intervals, had to be integrated with the large shelf-scale structural model. In order to achieve this with the limited and localized data, some globalization and stationarity assumptions were made throughout the geostatistical modeling process as described below. The resulting vertical range of 14 m for the exponential variogram model is consistent with the observed thickness of lithologic units in the borehole logs identified by (Mountain et al. 2020). This represents the smallest scale of vertical heterogeneity that will be captured by this model.
Standard variogram analysis employed in this study depends on two-point statistics, therefore, it is a measure of linear continuity and cannot efficiently capture curvilinear features (Hashemi et al. 2014). The clinothems imaged in the mid-Miocene sequences display a characteristic sigmoidal shape, most pronounced in sequences m54 and m58 (Fig. 4). By choosing a single rotational angle per direction for each sequence, we simplify the represented geometries. As a result, our analysis may underestimate the property correlation length. A review of possible frameworks for second-order non-stationary modeling to capture more complex geometries has been given by Fouedjio (2017). These include locally-varying anisotropies (Xu 1996;Martin et al. 2019), spatial transformation (Sampson and Guttorp 1992;Mallet 2002) and Non-Euclidean distance approaches (Boisvert et al. 2009). An integration of these concepts might improve the representation of shelf geometries in the final model. An alternative solution would be to calculate cross-well variograms (Deutsch and Pyrcz 2014), but such an approach would require more borehole data and a very detailed manual interpretation would be necessary for meaningful results. Considering the very low dip angle of the shelf platform, we assume that a linear approximation of the sigmoidal shapes is sufficiently representative.
We used a simple exponential variogram model with nugget component for all directions to avoid over-fitting the variogram to a small part of the full model domain. More sophisticated nested directional variogram structures might be able to better capture features like cyclicity in the vertical variogram (Gringarten and Deutsch 2001) (e.g. hole effect model), allowing improved modeling of small-scale variance. More data or additional expert knowledge would be required to define such a complex model of spatial correlation. A solid reasoning that this would be valid for the full shelf domain would be required, and a further spatial subdivision of the domain might be necessary, such as considering individual systems tracts.
In order to incorporate a horizontal trend in the model domain, we assumed that average of porosity values over the well log varies as a function of the well position in the model domain. We assume that the offshore fining trend, which is characteristic of shoreface sediments, will correlate with the average porosity trend. In a region-scale reservoir characterization study, Amthor and Moun (1994) used a cross-plot of well-scale averaged porosity versus geographical location to derive an empirical trend from 18 wells within a target formation. Average porosity values of wells across the entire study location also fell along this regression, suggesting that regional porosity prediction might be possible using this approach. We calculate a horizontal trend using only 3 wells as control points, resulting in a very gradual linear trend of increasing porosity offshore. Based on core analysis reports presented by Mountain et al. (2020), shale sediments showed a higher porosity range than sands. A gradual offshore increase in porosity is consistent with this, therefore, we consider this trend to be plausible, albeit, not very well constrained.
We propose that these assumptions and simplifications are reasonable to create a representative model that captures geological heterogeneity in 3D. While small-scale variability in subsurface properties might not be captured, we are able to represent the heterogeneity within the seven stratigraphic units, while respecting the governing geometries of the shelf setting, and the available porosity data as constraints. The resulting anisotropy ratios between the two horizontal-and the vertical directional variogram are within the expected range for deltaic environments (i.e. on the order of 100:1, (Deutsch and Pyrcz 2014)). The correlation lengths obtained for porosity in each sequence are consistently higher in the along-shore direction than the dip direction. This suggests elongated reservoir units, consistent with a wave-dominated depositional environment. It is important to characterize this correlation between the geometry of the depositional model and flow property distribution, as this is relevant to understanding flow and transport processes on a the continental shelf.
Variogram modeling is sensitive to several parameters, namely the maximum lag distance, bin width and the chosen theoretical model. For directional variograms the bandwidth and the tolerance angle affect the result. For this study, parameters were chosen based on literature suggestions (Deutsch and Pyrcz 2014;Gringarten and Deutsch 2001). For robust evaluation, a threshold for minimum number of point pairs was set based on the visual inspection of resulting experimental variograms. Further studies can benefit from a focus on sensitivity analysis to further constraint the suite of stochastic models. Note also that other geostatistical methods to integrate secondary seismic data like cokriging or collocated cokriging (Chiles and Delfiner 2012) were not considered in this work, as secondary data is only available along few seismic sections.
The Randomization method used in this work for the random field generation implemented in GSTools (Müller and Schüler 2021) has been successfully applied in hydrogeological studies (Comunian and Giudici 2021;Zech et al. 2021). The quality of the generated fields has been proven in theoretical comparisons with other methods (Kramer et al. 2007). Heße et al. (2014) highlighted the suitability of this method for reproducing small-scale variability within the generated field. The final model has a resolution of 500 m x 500 m x 20 m, corresponding to the resolution of the underlying geo-structural model. It has to be noted that it is straightforward to calculate higher-resolution realizations using the provided code if required for specific purposes.

Flow model quality and uncertainty
A quantitative assessment of the model realizations would require numerical flow simulation to compare the model performance to observed field data. This paper focuses on presenting the model paramterization workflow and making the data openly available. Here we make an initial qualitative analysis of the porosity and permeability field. This is based on geological knowledge of the depositional model. A basic preliminary quality check is presented in Fig.13. In this plot we can see that the simulated values match the overall trends of the observed porosity data, capturing the local variability of the porosity in the sediment column. The blind test realizations demonstrated that the model still matched the observed porosity at well site M28, which was removed from the conditioning data. The improvement in the fully conditioned realization shows that the shelf model can be further improved by the addition of data from future drilling campaigns.
The model is increasingly uncertain with distance away from the well locations, and in the absence of additional well data, these values cannot be cross-checked. However, the range of porosity values from 6 % -80 %, and the permeability values from 10 12 m 2 -10 18 m 2 , are entirely consistent with field observations to date. The porosity permeability relationship is constrained by ten laboratory measurements obtained mostly in sandy samples. As a result, more detailed modeling of permeability that accounts for grain size, such as the Kozeny-Carmen relation (Kozeny 1927), distribution could not be considered.

Conclusions
This study presents an integrated modeling approach with specific consideration for incorporating key components of the depositional environment into a flow property model. The combination of a seismic sequence stratigraphic interpretation with a comprehensive geostatistical analysis results in a geologically representative model that captures the heterogeneity of the shelf sediments in 3D, despite the limited well log data available. The resulting model highlights how anisotropy of flow properties are correlated with the geometries of the depositional model on the New Jersey shelf.
Using high-quality seismic profiles we were able to integrate a pseudo acoustic-impedance seismic attribute into a variogram modeling workflow through a pixel-based analysis. This study demonstrates that incorporating the seismic attribute in this way, is an effective and computationally inexpensive approach to overcome the challenge of sparsely distributed data, and to reliably estimate subsurface porosity distribution in the study area.
The modeling approach, integrates small-scale porosity data from core analysis to a shelf-scale model, while accounting for porosity compaction as well as basinward trends. The combination of a seismic sequence stratigraphic interpretation with a comprehensive geostatistical analysis results in a geologically representative model that captures the heterogeneity of the New Jersey shelf sediments in 3D. The current model does not capture the curvilinear geometries associated with stratigraphic clinoforms, instead, these are more simply represented as dipping linear trends. However, it is well suited for numerical investigation of the influence of geological heterogeneity on flow dynamics in the multi-layered offshore freshened groundwater system and other processes on the New Jersey continental shelf. The utilization of open-source tools (GemPy, GSTools) ensures that the model can be reproduced for future studies using supporting materials.

Page 17 of 19 11
in the preparation of this manuscript. We also thank the anonymous reviewers whose comments added to the value of this paper.
Funding Open Access funding enabled and organized by Projekt DEAL. Funding for this project is provided by the Deutsche Forschungsgemeinschaft (DFG) within the Priority Programme 527 -"International Ocean Discovery Program" (IODP) under Grants RE-3863/2-2 and BU1364/15-1.

Availability of data and materials
The data and code used to generate this model has been made available under a GNU General Public License v3.0 and can be accessed on the Zenodo open-access repository (https:// doi. org/ 10. 5281/ zenodo. 63274 50).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.