Thermal effects on flow and salinity distributions in heterogeneous coastal aquifers with fixed-flux inland boundaries

Freshwater-seawater (FW-SW) temperature contrasts widely exist in natural coastal aquifers. The significant effects of thermal forcing on water flow and salinity distributions for homogeneous aquifers have been demonstrated recently, however, the impact on heterogeneous aquifers remains unclear. This study conducted simulations of variable-density flow, and heat and salt transport with Monte-Carlo realizations of log-normally distributed permeability fields to examine such impacts. The averaged results showed that warmer freshwater could lead to a significant landward intrusion of freshwater-seawater interface in the heterogeneous aquifer. The random permeability fields increased the thermal effects of warmer freshwater and thus facilitated landward seawater intrusion. Furthermore, under warmer seawater conditions, salt dispersion was enhanced and density effects were reduced in heterogeneous coastal aquifers, thus altering the two opposing seawater circulation cells induced by double diffusion of salt and heat. The clockwise seawater circulation was strengthened whereas the anticlockwise one was weakened. Sensitivity analyses showed that an increased variance of permeability field further inhibited the anticlockwise seawater circulation cell caused mainly by heat diffusion. A larger correlation length of permeability field facilitated the thermal effect on the salinity distribution, increasing the associated uncertainty range caused by FW-SW temperature contrasts.

Key points:

➢ Freshwater-seawater interfaces are significantly affected by freshwater temperature in heterogeneous aquifers

➢ Salinity distributions are affected by variances and correlation lengths of permeability fields

➢ Double diffusion of salt and heat is inhibited with a larger variance of permeability field

Introduction

Flow, freshwater-seawater (FW-SW) mixing and associated biogeochemical reactions in coastal aquifers affect both marine and terrestrial ecosystems (Moore and Joye, 2021; Santos et al., 2021). With a terrestrial hydraulic gradient, inland fresh groundwater flows into the sea (Burnett et al., 2003). The FW-SW density contrast induces the convective circulation of seawater below lighter freshwater (Cooper, 1959; Smith, 2004). The two behaviors could form a FW-SW transition zone in the coastal aquifer, which could be treated as a sharp freshwater-seawater interface (FSI) when the transition zone is narrow relative to the aquifer thickness (Bear, 1979). As the freshwater recharge flux is large, density-dependent effects are weakened and the FSI is pushed seaward (Chang and Clement, 2012; Abdoulhalik and Ahmed, 2018).

FW-SW temperature contrasts have been widely reported around global coastlines (Wilson, 2005; Befus et al., 2013; Paldor et al., 2020; Pu et al., 2020). Temperature could alter the fluid density and viscosity in coastal aquifers, regulating density-driven flow and FW-SW mixing processes (van Lopik et al., 2015). Salinity distributions and groundwater discharge could thus be significantly regulated as the thermal regime changes in coastal aquifers (Nguyen et al., 2020; Pu et al., 2020). Wilson (2005) pointed out that geothermal convection could be developed in regional-scale coastal flow systems, contributing to the saline submarine groundwater discharge. The geothermal gradient modifies salinity distributions, particularly for the FSI toe (Blanco-Coronas et al., 2021). Nguyen et al. (2020) examined the impact of different seawater temperatures on the tidally influenced coastal unconfined aquifer. They found that higher seawater temperatures significantly increase tide-induced seawater circulations and thus affect fresh groundwater discharge zones and saltwater wedges. Recently, Pu et al. (2020) revisited the classical Henry problem for coastal confined aquifers and considered FW-SW temperature contrasts. Both laboratory experiments and numerical simulations demonstrated that warmer freshwater or colder seawater can lead to landward FSI in homogeneous coastal aquifers. Furthermore, due to double diffusion of salt and heat, a higher seawater temperature, once reaching certain levels, could cause two seawater circulation cells in the saltwater wedge (clockwise and anticlockwise circulations). The formations of two opposing circulation cells remarkably increase total fluxes of seawater circulations in homogeneous coastal aquifers (Pu et al., 2020).

While previous studies have demonstrated significant effects of thermal forcing on water flow and salinity distributions, geological heterogeneities were overlooked (Nguyen et al., 2020; Pu et al., 2020, 2021). The heterogeneity with spatial variability is the inherent property of coastal aquifer systems (Werner et al., 2013; Micallef et al., 2021), which may alter groundwater flow and strengthen the salt mixing and spreading (Kim et al., 2006; Michael et al., 2016). With the assumption of isothermal conditions, considerable efforts have been devoted to evaluating the effect of aquifer heterogeneity on salinity distributions and water circulations (e.g., Geng and Michael, 2020; Yu and Michael, 2019; Lu et al., 2013; Guo and Li, 2015). The random permeability field is a typical geostatistical method for simulating spatial heterogeneity, allowing a more realistic representation of complex geological properties (Freeze, 1975; Bellin and Rubin, 1996). Abarca et al. (2007) examined the Henry problem in heterogeneous aquifers with random permeability fields. The study demonstrated that salt dispersion and averaged permeability are predominant factors for the extent of seawater intrusion. Geological heterogeneity leads to increased salt dispersion around salt fronts and thus weakens density gradients, which causes the FSI to retreat seaward (Mahmoodzadeh and Karamouz, 2019; Ketabchi and Jahangir, 2021). Pool et al. (2015) examined the combined effects of heterogeneity and tides on the FSI and found that random heterogeneity enhanced the tide-induced salt mixing and spreading. Kreyns et al. (2020) conducted 3D density-dependent groundwater flow simulations to investigate the effects of heterogeneity on coastal volcanic aquifers and highlighted that the complex salinity distributions in heterogeneous aquifers enhance density-driven saltwater circulations increasing the ratio of saline to fresh groundwater discharge. Michael et al. (2016) demonstrated that combined effects of geological heterogeneity and density gradients could drive multiple circulation cells and result in more significant saltwater circulation rates. Yu and Michael (2022) also pointed out that the simplified geological model underestimates density-driven saltwater circulations. Geng et al. (2020a) found that in tidally influenced coastal beach aquifers, heterogeneity significantly alters groundwater flow in the intertidal circulation zone and leads to highly variable transit time. Although these studies revealed significant impacts of geological heterogeneity on water flow and salinity distributions in coastal aquifers, how geological heterogeneity behaves under varying thermal regimes remains unclear.

In this study, we conducted Monte-Carlo simulations considering random permeability fields and three FW-SW temperature contrasts to address the following research questions: 1) Whether salinity distributions and water circulations in heterogeneous coastal aquifers caused by temperature contrasts are consistent with those in homogeneous cases? 2) how do heterogeneous permeability fields change the thermal effects? 3) how do the variances and correlation lengths of permeability fields regulate salinity distributions considering temperature contrasts?

Conceptual model and simulations

The Henry problem (Henry, 1964) is a representative model of density-dependent flow and salt transport in natural coastal aquifers. Therefore, a similar conceptual model and setup have been adopted. As shown in Figure 1, a two-dimensional coastal confined aquifer (with a 400 m length and a 100 m height) perpendicular to the shoreline was considered. According to the setup of boundary conditions in the Henry problem (Henry, 1964), we adopted a fixed freshwater flux of 5.7 m³/m/d for the landward boundary CD. The seaward boundary AB was assigned to a fixed head of 100 m seawater level. Boundaries BC and AD are confined by impermeable layers, i.e., assumed to be no-flow boundaries (Henry, 1964). The seawater salinity was set to 35 ppt. Note that the conceptual model has been widely adopted in previous studies (e.g., Croucher and O’Sullivan, 1995; Simpson and Clement, 2004; Abarca et al., 2007; Voss and Provost, 2008; Sebben et al., 2015).

FIGURE 1

FIGURE 1. The conceptual model of a heterogeneous coastal confined aquifer with domain (ABCD) and boundary conditions indicated. The salinity and temperature are fixed at both the landward and seaward boundaries. The x-z coordinate origin is set at Point A (bottom right side of the confined aquifer).

The simulations were carried out with SUTRA-MS (Hughes and Sanford, 2005), a 2-D groundwater flow model coupled with solute and heat transport. In this study, the saturated groundwater flow in confined aquifers is based on Darcy’s law, as follows:

$\varphi \left(\frac{\partial {\rho }_f}{\partial C}\frac{\partial C}{\partial t}+\frac{\partial {\rho }_f}{\partial T}\frac{\partial T}{\partial t}\right)=-\nabla \cdot \left({\rho }_f\overrightarrow{q}\right)+{\rho }_F{Q}_{\mathrm{F}}(1)$

$\overrightarrow{q}=-\left(\frac{k}{\mu }\right)\cdot \left(\nabla p-{\rho }_{f}g\right)(2)$

where $\varphi$ is the soil porosity [−]. $C$ and $T$ are the fluid salinity [-] and temperature [°C], respectively. $t$ represents the time [T] and $\overrightarrow{q}$ represents the Darcy velocity [LT⁻¹]. $Q_{\mathrm{F}}$ and ${\rho }_F$ are the fluid source [T⁻¹] and the corresponding density of fluid source [ML⁻³], respectively. $\mu$ and $g$ are the fluid viscosity [ML⁻¹T⁻¹] and the gravitational acceleration [LT⁻²], respectively. $k$ and $p$ represent the intrinsic permeability [L²] and the pore-water pressure [ML⁻¹T⁻²], respectively. ${\rho }_f$ is the fluid density [ML⁻³]. Note that the fluid density is assumed to be linearly associated with the temperature and salinity as done in Pu et al. (2020), i.e., ${\rho }_f=1002+0.7705\mathrm{C}-0.2205\mathrm{T}$. Fluid viscosity is negatively related to fluid temperature ($\mu =239.4\times {10}^{\left(\frac{248.37}{T+133.15}\right)-7}$, given by (Hughes and Sanford, 2005)).

Salt and heat transports are described by advection-dispersion equations, as follows:

$\frac{\partial \left(\varphi {\rho }_{f}C\right)}{\partial t}=-\nabla \cdot \left({\rho }_f\overrightarrow{q}C\right)+\nabla \cdot \left(\varphi {\rho }_f{\mathbf{D}}_S\nabla C\right)+{\rho }_F{Q}_F{C}_F(3)$

$\frac{\partial \left[\varphi {\rho }_f{c}_f+\left(1-\varphi \right){\rho }_s{c}_s\right]T}{\partial t}=-\nabla \cdot \left({\rho }_f{c}_f\overrightarrow{q}T\right)+\nabla \cdot \left({\rho }_f{c}_f{D}_T\nabla T\right)+c_f{\rho }_F{Q}_F{T}_F(4)$

$D_T=\frac{\varphi {\lambda }_f+\left(1-\varphi \right){\lambda }_s}{{\rho }_f{c}_f}(5)$

where ${\mathbf{D}}_S$ is the diffusion coefficient tensor (Voss and Provost, 2008). $C_F$ and $T_F$ represent the salinity [-] and temperature [°C] of fluid source, respectively. $c_s$ and $c_f$ are the solid matrix specific heat [EM⁻¹°C⁻¹] and the fluid specific heat [EM⁻¹°C⁻¹], respectively. ${\rho }_s$ is the solid matrix density [ML⁻³]. ${\lambda }_s$ and ${\lambda }_f$ are the thermal conductivity [ET⁻¹L⁻¹°C⁻¹] of solid matrix and fluid, respectively. $D_T$ is the effective thermal diffusivity [L²T⁻¹]. In Eq 4, 5, $\varphi {\rho }_f{c}_f+\left(1-\varphi \right){\rho }_s{c}_s$ and $\varphi {\lambda }_f+\left(1-\varphi \right){\lambda }_s$ represent the bulk heat capacity and bulk thermal conductivity of the saturated soil, respectively.

A log-normal permeability field was chosen to mimic geological heterogeneity, as widely adopted previously (e.g., Freeze, 1975; Abarca et al., 2007; Geng et al., 2020b; Ketabchi and Jahangir, 2021). The random permeability field is generated using Hydro_gen (Bellin and Rubin, 1996). For the base case, the permeability distribution was characterized by a mean logarithm value of −25.82 (corresponding to geometric mean permeability of 10⁻¹¹ m²) and a logarithm variance (${\sigma }_{\ln k}^2$) of 1. An exponential-decay covariance function is adopted to represent the spatially correlated structure of permeability field. The key parameters of random fields are the variance (${\sigma }_{\ln k}^2$) and directional correlation lengths (horizontal correlation length ${\lambda }_x$ of 0.5 m and vertical correlation length ${\lambda }_z$ of 0.05 m). In addition to the cases with randomly-distributed permeability, the homogeneous cases with uniform permeability of 10⁻¹¹ m² (similar to that used in Pu et al. (2020)) were also examined.

In this study, we investigate the effect of three FW-SW temperature contrasts on heterogeneous coastal aquifers ($T_s$ = 35°C and $T_f$ = 20°C, $T_s$ = 20°C and $T_f$ = 20°C, and $T_s$ = 20°C and $T_f$ = 35°C, Table 1). Note that there exists a temperature contrast larger than 15°C in the natural coastal system (Kohout, 1965; Hughes et al., 2007). The temperature contrast of 15°C was also adopted in previous studies (Nguyen et al., 2020; Pu et al., 2020, 2021). Furthermore, we conducted sensitivity analyses on log-permeability variances (${\sigma }_{\ln k}^2$ = 1, 2, and 3, Table 2 and Supplementary Figure S1) and field correlation lengths (${\lambda }_x$ = 0.5, 2.5, and 12.5 m, Table 3 and Supplementary Figure S2). The correlation length ratio of z-direction to x-direction was set to 1/10 in our studies (Gelhar, 1986).

TABLE 1

TABLE 1. Key results for homogeneous and heterogeneous cases influenced by different FW-SW temperature contrasts.

TABLE 2

TABLE 2. Key results for heterogeneous cases with different permeability variances and FW-SW temperature contrasts.

TABLE 3

TABLE 3. Key results for heterogeneous cases with different correlation lengths and FW-SW temperature contrasts.

Other parameter values used in this study are consistent with those adopted in the large-scale model considered by Pu et al. (2020). The porosity was set to 0.35. The longitudinal dispersivity (${\alpha }_L$) and transverse dispersivity (${\alpha }_T$) were set to 1 and 0.1 m, respectively (Hunt et al., 2011). The molecular diffusivity was assigned to 10^–9 m²/s (Voss and Provost, 2008). The fluid specific heat and solid matrix specific heat were set to 4,182 and 840 J/kg/°C, respectively (Engineering ToolBox, 2003a). The fluid thermal conductivity and solid thermal conductivity were set to 0.6 and 3.5 J/m/°C/s, respectively (Engineering ToolBox, 2003b; Anderson, 2005).

A model mesh was generated with a uniform discretization of 1 m in both x-direction and z-direction. This spatial discretization ensures numerical stability requirements (i.e., Peclet criteria), as suggested by Voss and Provost (2008). 50 Monte-Carlo simulations of randomly-distributed permeability were conducted for each model setup, leading to 750 simulation runs. All simulations were run to the steady state, i.e., pressure, salinity, and temperature unchanged.

Results

Model validation at the laboratory scale

Firstly, we validated the numerical model against experimental results (Pu et al., 2020) under FW-SW temperature contrasts ($T_s$ = 40°C and $T_f$ = 25°C, $T_s$ = 25°C and $T_f$ = 25°C, $T_s$ = 25°C and $T_f$ = 33°C) and heterogeneous permeability fields. The model parameters and boundary conditions of homogeneous aquifers were similar to those adopted in Pu et al. (2020). For heterogeneous confined aquifers, the geometric mean permeability of 5 × 10⁻¹⁰ m² and the logarithm variance of 0.5 were used (the randomly-distributed permeability is shown in Supplementary Figure S3).

As expected, the randomly-distributed permeability impacted the extent of FW-SW mixing, changing groundwater flow paths and steady-state salinity distributions (Figure 2). However, the associated changes in salinity distributions due to random permeability fields were not significant at the laboratory scale. The simulated steady-state FSI (i.e., 50% seawater isohalines of 17.5 ppt) for heterogeneous aquifers was also in well agreement with the measured results (Pu et al., 2020). Differences of FSI toe locations (intersection of the 50% seawater isohaline with the aquifer basement) between heterogeneous simulations and flume experiments were 0.13, 0.17, and 0.14 m, respectively, for cases with warmer seawater, isothermal condition and warmer freshwater (Figures 2B,D,F). This demonstrates the effectiveness of simulations in heterogeneous coastal aquifers subjected to temperature contrasts and could be used as model validation.

FIGURE 2

FIGURE 2. The laboratory-scale salinity distributions under different FW-SW temperature contrasts. The left-hand (A, C, E) and right-hand (B, D, F) panels are the results for homogeneous and heterogeneous aquifers, respectively. The black solid and dashed lines are, respectively, for the 50% seawater concentration isohalines (17.5 ppt) and the measured (Pu et al., 2020) freshwater-seawater interface (FSI). The white lines show flow streamlines. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

Thermal effects on salinity distributions and water circulations

Salinity distributions of heterogeneous aquifers under different FW-SW temperature contrasts ($T_s$ = 35°C and $T_f$ = 20°C, $T_s$ = 20°C and $T_f$ = 20°C, $T_s$ = 20°C and $T_f$ = 35°C) were examined here. For the isothermal condition, the location of steady-state FSI toe in a Monte-Carlo simulation was −113 m, different from that for the homogeneous aquifer (x = −174 m, Supplementary Figure S4). Furthermore, the steady-state FSI and groundwater flow paths were unsmooth in the heterogeneous case. To further quantify the changes in heterogeneous aquifers, we calculated the averaged results (e.g., salinity distribution and water flux) of 50 heterogeneous permeability simulations. The averaged results smoothed out the disturbance of randomly-distributed permeability to salinity distributions and flow paths (Figure 3). Therefore, the variations caused by thermal forcing could be exhibited clearly for heterogeneous coastal aquifers.

FIGURE 3

FIGURE 3. Salinity distributions in coastal aquifers subjected to different FW-SW temperature contrasts. The left-hand (A, C, E) and right-hand (B, D, F) panels are the results for homogeneous and heterogeneous aquifers, respectively. The results for heterogeneous aquifers (${\sigma }_{\ln k}^2$ = 1 and ${\lambda }_x$ = 0.5 m) are averaged (n = 50). The black dashed lines are the 50% seawater concentration isohalines (17.5 ppt), i.e., the FSI. The white lines show flow streamlines. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

With warmer seawater, the averaged FSI toe slightly retreated seaward (from x = −112 to −102 m) for heterogeneous coastal aquifers, consistent with the results from homogeneous cases (Figure 3). The total amount of salt stored in the aquifer ($SM$) accordingly decreased (from 38,693 to 32,559 kg/m, Table 1). Similar to homogeneous aquifers, warmer seawater decreased seawater density (from 1024.56 to 1021.25 kg/m³, Figures 4B,D), and thus weakened pressure gradients along the FSI in coastal confined aquifers (Glover, 1959). With warmer seawater, groundwater heads at landward boundaries were slightly smaller than that under the isothermal conditions for both homogeneous and heterogeneous aquifers (103.89 compared with 104.05 m and 105.20 compared with 105.37 m, respectively). This caused a slight decrease in the lateral hydraulic gradient. On the whole, the effect of density gradient on the FSI was more significant than that of the lateral hydraulic gradient, resulting in the FSI retreating seaward.

FIGURE 4

FIGURE 4. Density distributions in coastal aquifers subjected to different FW-SW temperature contrasts. The left-hand (A, C, E) and right-hand (B, D, F) panels are the results for homogeneous and heterogeneous aquifers, respectively. The results for heterogeneous aquifers (${\sigma }_{\ln k}^2$ = 1 and ${\lambda }_x$ = 0.5 m) are averaged (n = 50). The black dashed lines are the 50% seawater concentration isohalines (17.5 ppt), i.e., the FSI. Note that only the sections from x = −300 to 0 m are presented. The white lines show flow streamlines. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

In contrast, warmer freshwater facilitated saltwater intrusion into the heterogeneous freshwater aquifer (Figures 3D,F). The averaged FSI moved landward with the toe changing from x = −112 to −170 m. Correspondingly, the associated $SM$ increased remarkably (from 38,693 to 67,340 kg/m). Warmer freshwater decreased the freshwater density and thus increased the FW-SW density difference (with an increase of 3.31 kg/m³). This favored the landward intrusion of FSI compared to the isothermal case. More importantly, the higher freshwater temperature can largely decrease the water viscosity in the freshwater zone and thus increase the hydraulic conductivity (approximately increased by 50%), declining the lateral hydraulic gradient for the flux-controlled system (Pu et al., 2020). The lateral hydraulic gradient (calculated by the ratio of total freshwater head difference between landward and seaward aquifer boundaries to the horizontal aquifer length) decreased from 0.0134 to 0.0111 for the heterogeneous aquifer with warmer freshwater. Therefore, warmer freshwater caused the saltwater wedge to advance landward remarkably. This was consistent with the thermal impacts on homogeneous aquifers (Pu et al., 2020).

Compared with homogeneous aquifers, the averaged FSI in heterogeneous aquifers was closer to the seaward boundary (with smaller $SM$, Table 1). The retreated FSI was attributed to enhanced lateral hydraulic gradients in heterogeneous aquifers. The averaged lateral hydraulic gradients of heterogeneous aquifers were larger than that of homogeneous cases for the three temperature contrasts (0.0130 compared with 0.0097, 0.0134 compared with 0.0101, and 0.0111 compared with 0.0087, respectively). Similarly, randomly-distributed permeability also led to a relatively narrow temperature mixing zone for both warmer seawater and freshwater cases (Supplementary Figure S5). Furthermore, warmer freshwater significantly contributed to the hydraulic conductivity of heterogeneous aquifers. The groundwater head at the landward boundary decreased by 0.95 m, larger than that for the homogeneous aquifer (0.58 m, Supplemntary Figure S6). Warmer freshwater largely decreased the lateral hydraulic gradient compared with the homogeneous aquifer (difference of 0.0023 compared with 0.0014). The averaged relative change of $SM$ caused by warmer freshwater for the heterogeneous aquifer was more significant than that for the homogeneous aquifer (74% compared with 60%). The results suggest that randomly-distributed permeability enhanced the impacts of warmer freshwater on hydraulic gradients and salinity distributions.

For the isothermal case, mechanical dispersion and molecular diffusion could cause the salt to move towards the freshwater aquifer. That induced a FW-SW mixing zone and contributed to density-driven convective circulations in the saltwater wedge (Smith, 2004). With randomly-distributed permeability, salt dispersion in the saltwater wedge was enhanced (Kerrou and Renard, 2010). It caused a larger seawater circulation flux (0.47 m³/m/d). The relative change caused by random permeability field was 67% for the total seawater influx, compared with uniform permeability. Further, random permeability fields changed water fluxes per unit area across seaward boundaries. Water exchange across the aquifer-sea interface was larger for high-permeability zones of heterogeneous aquifers. And the freshwater efflux curves were no longer smooth and monotonous (Figure 5B). The maximum freshwater efflux per unit area for the heterogeneous aquifer was much less than that for the homogeneous aquifer. It could be attributed to dissipated advective forces as permeability randomly distributed in coastal aquifers.

FIGURE 5

FIGURE 5. Per-unit-area water flux across seaward boundary in (A) homogeneous and (B) heterogeneous confined aquifers. The results of different FW-SW temperature contrasts are indicated by lines with different colors. The results for heterogeneous aquifers (${\sigma }_{\ln k}^2$ = 1 and ${\lambda }_x$ = 0.5 m) are averaged (n = 50). $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

With warmer seawater, we also found the two opposing circulation cells in saltwater wedges of heterogeneous aquifers (Figure 3B). FW-SW temperature contrasts increased the complexity of density gradients and regulated density-driven seawater circulations. As the effective thermal diffusivity (5.93 × 10⁻⁷ m²/s) was significantly larger than the salt diffusivity (10⁻⁹ m²/s), double diffusion of heat and salt within the saltwater wedge induced the formation of two seawater circulation cells in opposing directions. It was consistent with the homogeneous results (Pu et al., 2020). Under warmer seawater conditions, flow velocities in the upper (clockwise) circulation cells were largely increased, influencing total seawater influxes. Although the clockwise seawater circulation zone was significantly smaller than that for the isothermal case (x = 0 and z = 30–67 m in comparison with x = 0 and z = 0–70 m), the averaged amount of seawater influx for the clockwise circulation cell was close to the total seawater influx of isothermal case (0.47 m³/m/d, Table 1). The anticlockwise seawater influx induced mainly by heat diffusion was 0.07 m³/m/d. The average percentage of anticlockwise circulating seawater influx to total seawater influx (hereinafter termed the percentage of reverse seawater influx) was 15%, less than half of that for the homogeneous aquifer (34%). Due to intensified salt dispersion in heterogeneous aquifers, the upward convection was strengthened whereas the downward convection mainly induced by heat diffusion was weakened, i.e., randomly-distributed permeability weakened the effects of buoyancy gradients induced by heat diffusion.

With warmer freshwater, larger FW-SW density differences enhanced upward buoyancy forces of saltwater around salt fronts, increasing density-driven circulation rates of heterogeneous aquifers. Further, warmer freshwater heated part of saltwater within the saltwater wedge through heat conduction, thus influencing the density distribution (Figure 4E). This caused saltwater flow streamlines to tilt upwards further, regulating the circulation paths. The seawater circulation zone significantly expanded for the heterogeneous aquifer with warmer freshwater (x = 0 and z = 0–80 m, x = 0 and z = 0–70 m, respectively, for warmer freshwater and isothermal cases). It produced a total seawater influx of 0.61 m³/m/d (30% relative increase) for the warmer freshwater case. Compared with the homogeneous aquifer, the effect of warmer freshwater on the variation of total seawater influx was also weakened in the heterogeneous aquifer (relative change, 30% compared with 61%).

In addition, we applied box plots to show the distributions of quantitative data for toe locations, $SM$ and seawater influxes of 50 heterogeneous simulations (Figure 6). For isothermal cases, the toe locations varied between −115 and −107 m and the $SM$ ranged between 36,522 and 41,106 kg/m. Thermal effects could influence uncertainty ranges caused by randomly-distributed permeability. For toe locations and $SM$, the uncertainty range was narrowed by warmer seawater and enlarged by warmer freshwater. For example, warmer seawater led to a decreased interquartile range for the $SM$, in comparison with the isothermal case (1,175 in comparison with 1,723 kg/m). In contrast, warmer freshwater further increased the interquartile range (1,723 in comparison with 1,850 kg/m). That could be attributed to the changed lateral hydraulic gradients and density differences influenced by thermal effects.

FIGURE 6

FIGURE 6. Box plots (n = 50) for the heterogeneous aquifers subjected to different FW-SW temperature contrasts. The left-hand, middle and right-hand panels are the results for (A) the toe location, (B) the total amount of salt mass stored in the aquifer (SM), and (C) the total seawater influx across the aquifer-sea interface, respectively. The dashed lines are the results for homogeneous aquifers. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

The uncertainty range of total seawater influx in heterogeneous aquifers was relatively large (Figure 6C). The minimum seawater influx was lower than that for the homogeneous aquifer (0.23 compared with 0.28 m³/m/d). In contrast, the maximum seawater influx of heterogeneous aquifer reached 1.05 m³/m/d, about four times that of homogeneous aquifer. This suggests that the random permeability field induced a significant variation of seawater circulation flux. As temperatures of seawater or freshwater increased, the uncertainty of total seawater influx was changed (Figure 6C). With warmer seawater, the lower quantile of total seawater influxes was overlapped with the value of homogeneous aquifer. The percentages of reverse seawater influxes varied from 5.53% to 26.84% for heterogeneous cases, suggesting that random permeability fields could change two opposing seawater circulation cells significantly. With warmer freshwater, the range (maximum-minimum) of total seawater influx was larger than that of the isothermal case (1.03 compared with 0.81 m³/m/d). Whereas the interquartile range of warmer freshwater was smaller than the isothermal case (0.18 compared with 0.24 m³/m/d). The results demonstrate the importance of thermal effects on heterogeneous aquifers and suggest that permeability fields could significantly regulate the impacts.

Sensitivity analysis on variances of permeability fields

In this section, we conducted additional simulations on variances of permeability fields (${\sigma }_{\ln k}^2$ = 1, 2, and 3, respectively) to further examine the variations of heterogeneous aquifers caused by FW-SW temperature contrasts (n = 50 for each setup). The correlation length and other parameters remained unchanged in these numerical simulations. The increased ${\sigma }_{\ln k}^2$ led to an extended range of randomly-distributed permeability (Supplementary Figure S1), influencing groundwater flow and salinity distributions. As ${\sigma }_{\ln k}^2$ increased, the impact of geological heterogeneity was significantly enhanced. This resulted in a higher lateral hydraulic gradient (0.0134, 0.0188, and 0.0274, respectively, for ${\sigma }_{\ln k}^2$ = 1, 2, and 3 in the isothermal case). Therefore, the intrusion of FSI was remarkably inhibited for larger ${\sigma }_{\ln k}^2$ (for isothermal cases, the averaged toe locations moved from x = −112 to −66 and −34 m as ${\sigma }_{\ln k}^2$ increased from 1 to 2, and 3, Supplementary Figure S7). As expected, the thermal effects on salinity distributions were also significant for ${\sigma }_{\ln k}^2$ = 2 and 3. The larger ${\sigma }_{\ln k}^2$ increased relative changes caused by warmer seawater for the averaged $SM$ (relative decreases of 16%, 21% and 31%, respectively, for ${\sigma }_{\ln k}^2$ = 1, 2, and 3). Similarly, as ${\sigma }_{\ln k}^2$ increased, the variations caused by warmer freshwater were remarkably increased (relative increases of 74%, 108% and 186%, respectively, for the averaged $SM$ with ${\sigma }_{\ln k}^2$ = 1, 2, and 3). With the larger ${\sigma }_{\ln k}^2$, the decrease of lateral hydraulic gradient was more pronounced caused by warmer freshwater (from 0.0023 to 0.0039 and 0.0063). This led to changes in salinity distributions (Goswami and Clement, 2007).

The box plots were shown to compare thermal effects of heterogeneous aquifers with different permeability variances (Figure 8). The larger ${\sigma }_{\ln k}^2$ widened uncertainty ranges of toe locations due to the increased dispersion (Figure 8A). For isothermal conditions, the interquartile ranges of toe locations increased from 2.98 to 3.37 and 3.49 m as ${\sigma }_{\ln k}^2$ increased from 1 to 2 and 3, respectively. The uncertainty range of $SM$ tended to decrease as ${\sigma }_{\ln k}^2$ increased (Figure 8B). This was because the FSI significantly retreated seaward with the larger ${\sigma }_{\ln k}^2$. However, the associated differences in uncertainty ranges caused by different ${\sigma }_{\ln k}^2$ were not as evident as the averaged results for both the toe location and $SM$. The changes in uncertainty ranges caused by warmer freshwater were relatively larger than that caused by warmer seawater (Figures 8A,B), showing that warmer freshwater could significantly strengthen the variation of hydraulic conductivity and thus enlarge the differences in salinity distributions among heterogeneous simulations. The variations can be further enhanced by larger ${\sigma }_{\ln k}^2$. As ${\sigma }_{\ln k}^2$ increased from 1 to 2 and 3, the ratios of interquartile ranges (warmer freshwater to the isothermal condition) for $SM$ changed from 107% to 167% and 201%, respectively.

Concurrently, groundwater flow paths were adjusted by the variances of permeability fields (Figure 7). With the larger ${\sigma }_{\ln k}^2$, disturbances caused by aquifer heterogeneity to the flow and salt transport were intensified. The larger ${\sigma }_{\ln k}^2$ increased salt dispersion altering saltwater circulation rates and seawater influxes significantly. Although saltwater wedges retreated seaward as ${\sigma }_{\ln k}^2$ enhanced, the averaged seawater influxes remarkably increased (from 0.47 to 0.85 and 1.28 m³/m/d for isothermal cases with ${\sigma }_{\ln k}^2$ = 1, 2, and 3, Table 2). This suggests that local seawater circulation fluxes in freshwater discharge zones were significantly enhanced, resulting from geological heterogeneity and local density gradients (Supplementary Figure S8). Both warmer seawater and freshwater can increase total seawater influxes of heterogeneous aquifers with larger permeability field variances (Table 2). In particular, for warmer seawater, the anticlockwise circulation cell was much smaller as the ${\sigma }_{\ln k}^2$ increased. With ${\sigma }_{\ln k}^2$ = 3, no apparent anticlockwise seawater circulation cell demonstrated (Figure 7C). The associated averaged anticlockwise seawater influx was very small because salt dispersion was intensified and density effects (induced by temperature contrast) was weakened in highly heterogeneous aquifers. However, for ${\sigma }_{\ln k}^2$ = 3, the clockwise seawater influx under warmer seawater conditions reached 107% of the total seawater influx under isothermal conditions on average. This was due to enhanced seawater circulation rates for warmer seawater cases.

FIGURE 7

FIGURE 7. Averaged salinity distributions in heterogeneous aquifers under different permeability variances and FW-SW temperature contrasts (n = 50). The left-hand (A, D, G), middle (B, E, H) and right-hand (C, F, I) panels are the results for different variances of random permeability fields (${\sigma }_{\ln k}^2$ = 1, 2, and 3, respectively). ${\lambda }_x$ is fixed at 0.5 m. The black dashed lines are the 50% seawater concentration isohalines (17.5 ppt), i.e., the FSI. The white lines show flow streamlines. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

With larger ${\sigma }_{\ln k}^2$, the maximum total seawater influxes increased significantly for isothermal cases (1.05, 1.93, and 2.77 m³/m/d, respectively, for ${\sigma }_{\ln k}^2$ = 1, 2, and 3). The interquartile ranges of total seawater influxes increased from 0.24 to 0.36 and 0.66 m³/m/d as ${\sigma }_{\ln k}^2$ increased from 1 to 2 and 3, respectively. It indicates that the larger ${\sigma }_{\ln k}^2$ increased the uncertainty of total seawater influx. For ${\sigma }_{\ln k}^2$ = 2 and 3, both warmer seawater and freshwater could increase uncertainty ranges of total seawater influxes (Figure 8C). As stated above, the higher seawater and freshwater temperatures in heterogeneous aquifers could lower lateral hydraulic gradients, which favored salt dispersion induced by randomly-distributed permeability thus strengthening uncertainty of total seawater influxes. For warmer seawater, two circulation cells in the saltwater wedge were remarkably changed as ${\sigma }_{\ln k}^2$ increased. The percentages of reverse seawater influxes were lowered (Figure 9). The average percentages decreased from 15% to 6% and 2% as ${\sigma }_{\ln k}^2$ increased from 1 to 2 and 3, respectively. With ${\sigma }_{\ln k}^2$ = 3, most percentages of reverse seawater influxes were not significant in heterogeneous simulations (the upper quantile of 2%, n = 50). The results suggest that two opposing seawater circulation cells induced by double diffusion of salt and heat were very sensitive to and negatively correlated with the variance of heterogeneous permeability field.

FIGURE 8

FIGURE 8. Box plots (n = 50) for heterogeneous aquifers under different permeability variances (${\sigma }_{\ln k}^2$ = 1, 2, and 3, respectively) and FW-SW temperature contrasts. The left-hand, middle and right-hand panels are the results for (A) the toe location, (B) the total amount of salt mass stored in the aquifer (SM), and (C) the total seawater influx across the aquifer-sea interface, respectively. ${\lambda }_x$ is fixed at 0.5 m. The dashed lines are the results for homogeneous aquifers. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

FIGURE 9

FIGURE 9. Box plots (n = 50) for percentages of anticlockwise circulating seawater influx to the total seawater influx in heterogeneous aquifers with different permeability field variances (${\sigma }_{\ln k}^2$ = 1, 2, and 3, respectively). ${\lambda }_x$ is fixed at 0.5 m. The black dashed line is the result for the homogeneous aquifer.

Sensitivity analysis on correlation lengths of permeability fields

With other parameters unchanged, we tested the importance of correlation lengths (${\lambda }_x$ = 0.5, 2.5, and 12.5, respectively and ${\lambda }_z/{\lambda }_x$ = 0.1) for heterogeneous aquifers with random permeability fields (n = 50). The increased correlation length can reduce the difference among neighbored points and thus result in the formation of patches in the permeability field (Supplementary Figure S2). The larger high-permeability patch would produce preferential flow paths and modulate the flow and local FW-SW mixing in heterogeneous aquifers (Pool et al., 2015). An increase in ${\lambda }_x$ induced a slight decrease in the averaged lateral hydraulic gradient (0.0134, 0.0123, and 0.0115, respectively, for ${\lambda }_x$ = 0.5, 2.5, and 12.5 m in isothermal cases). Therefore, the FSI intruded landward slightly (Figure 10; Supplementary Figure S9). The averaged toe locations changed from x = −112 to −126 and −139 m as ${\lambda }_x$ increased from 0.5 to 2.5 and 12.5 m, respectively. With the increase of ${\lambda }_x$, the heterogeneous aquifer is likely to stratify densely with multiple patches. Thus, the wider FW-SW mixing zone for the averaged salinity distribution was found. Temperature contrasts further regulated the averaged salinity distributions in aquifers. Similarly, warmer freshwater caused a more significant variation than warmer seawater for ${\lambda }_x$ = 2.5 and 12.5 m, in comparison with isothermal cases (e.g., 74% compared with 18%, 73% compared with 17%, respectively, for the averaged $SM$ with ${\lambda }_x$ = 2.5 and 12.5 m).

Furthermore, we found that uncertainty ranges of toe location and $SM$ were significantly increased for heterogeneous aquifers with larger correlation lengths (Figures 11A,B). For isothermal conditions, as ${\lambda }_x$ increased from 0.5 to 2.5 and 12.5 m, the interquartile ranges of toe locations increased from 2.98 to 4.83 and 13.98 m, respectively. Correspondingly, the interquartile ranges of $SM$ increased from 1,723 to 3,243 and 10,373 kg/m. The variations of uncertainty ranges were attributed to the enhanced FW-SW mixing in heterogeneous aquifers. The uncertainty ranges of toe location and $SM$ were further widened as freshwater temperature increased (Figures 11A,B). With ${\lambda }_x$ = 12.5 m, the ratio of interquartile range (warmer freshwater to the isothermal condition) reached 117% for the toe location and 108% for the $SM$. The results demonstrate that the larger ${\lambda }_x$ facilitated effects of thermal forcing, thus increasing uncertainty ranges of salinity distributions in heterogeneous aquifers.

The impacts of different correlation lengths on salt dispersion were slight due to the similar permeability range of field. Thus, by averaging, the total seawater influx was similar (0.47, 0.46 and 0.47 m³/m/d, respectively, for ${\lambda }_x$ = 0.5, 2.5 and 12.5 m in isothermal conditions). The thermal effects on water circulations for ${\lambda }_x$ = 2.5 and 12.5 m were close to that of ${\lambda }_x$ = 0.5 m (Figure 10; Supplementary Figure S10). Both warmer seawater and freshwater can enhance the averaged seawater influxes in the cases with ${\lambda }_x$ = 2.5 and 12.5 m (Table 3). For example, with ${\lambda }_x$ = 12.5 m, warmer seawater and freshwater, respectively, led to relative increases of 9% and 38% for the total seawater influx. For warmer seawater, the different ${\lambda }_x$ slightly changed two opposing circulation cells induced by double diffusion of heat and salt (Figures 10A–C). The average percentages of reverse seawater influxes were 15%, 11%, and 12%, respectively, for ${\lambda }_x$ = 0.5, 2.5 and 12.5 m.

FIGURE 10

FIGURE 10. Averaged salinity distributions in heterogeneous aquifers under different correlation lengths and FW-SW temperature contrasts (n = 50). The left-hand (A, D, G), middle (B, E, H) and right-hand (C, F, I) panels are the results for different correlation lengths of permeability fields (${\lambda }_x$ = 0.5, 2.5, and 12.5 m, respectively). ${\sigma }_{\ln k}^2$ is fixed at 1. The black dashed lines are the 50% seawater concentration isohalines (17.5 ppt), i.e., the FSI. The white lines show flow streamlines. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

The larger ${\lambda }_x$ caused total seawater influxes to centralize slightly, narrowing uncertainty ranges (Figure 11C). The associated uncertainty range was slightly changed by seawater and freshwater temperatures for ${\lambda }_x$ = 2.5 and 12.5 m. With warmer seawater, the two circulation cells were affected and associated uncertainty ranges behaved differently for ${\lambda }_x$ = 2.5 and 12.5 m (Figure 12). The larger ${\lambda }_x$ was expected to reduce the uncertainty range for the percentage of reverse seawater influx. However, with ${\lambda }_x$ = 12.5 m, the uncertainty range of anticlockwise seawater circulation flux was increased, causing a larger difference in percentages of reverse seawater influxes (from 1% to 31%). The results showed that all percentages of reverse seawater influxes in heterogeneous aquifers were lower than that for the homogeneous aquifer (Figure 12).

FIGURE 11

FIGURE 11. Box plots (n = 50) for heterogeneous aquifers under different correlation lengths (${\lambda }_x$ = 0.5, 2.5, and 12.5 m, respectively) and FW-SW temperature contrasts. The left-hand, middle and right-hand panels are the results for (A) the toe location, (B) the total amount of salt mass stored in the aquifer (SM), and (C) the total seawater influx across the aquifer-sea interface, respectively. ${\sigma }_{\ln k}^2$ is fixed at 1. The dashed lines are the results for homogeneous aquifers. $T_s$ and $T_f$ are seawater and freshwater temperatures, respectively.

FIGURE 12

FIGURE 12. Box plots (n = 50) for percentages of anticlockwise circulating seawater influx to the total seawater influx in heterogeneous aquifers with different correlation lengths (${\lambda }_x$ = 0.5, 2.5 and 12.5 m, respectively). ${\sigma }_{\ln k}^2$ is fixed at 1. The black dashed line is the result for the homogeneous aquifer.

Discussion

The results on heterogeneous coastal confined aquifers showed that the freshwater-seawater interface intrudes further landward as freshwater temperature increases. This is consistent with the previous findings of homogeneous aquifers (Pu et al., 2020). For aquifers with randomly-distributed permeability, salinity distributions could be changed significantly under the same boundary conditions (Pool et al., 2015; Mahmoodzadeh & Karamouz, 2019). Yu and Michael (2022) applied the upscaling hydraulic conductivity methods and found that salinity distributions considering upscaling effects are completely different from the original one. However, our results suggested that the indicators of seawater intrusion varied within a relatively small range among Monte-Carlo simulations, compared with the mean value. For salinity distributions, differences in the averaged results caused by thermal forcing are easier to compare and generalize. Therefore, equivalent homogeneous model is still a feasible solution for investigating variations of salinity distributions in coastal aquifers affected by various forcing factors. For site-specific coastal aquifers, variations of salinity distributions are affected by combined specific geological architectures and hydrological forcing factors. Groundwater simulations considering spatial heterogeneity could be better for this situation.

Pu et al. (2020) pointed out that FW-SW temperature contrasts regulate salinity distributions by changing density differences and lateral hydraulic gradients in homogeneous coastal aquifers. For heterogeneous aquifers, permeability contrasts are critical factors influencing salinity distributions (Lu et al., 2013; Yu and Michael, 2022). Thermal effects on high-permeability zones are more significant than that on low-permeability zones. That is, the differences in salinity distributions caused by thermal effects would be enlarged for aquifers with bigger permeability contrasts (e.g., larger permeability field variances), particularly for impacts of freshwater temperatures. As freshwater temperature increases, the lateral hydraulic gradient of heterogeneous aquifer is further lowered, affecting the salinity distribution. In natural aquifer systems, geological heterogeneity with spatial variability is common. Due to heterogeneity, the effects of warmer freshwater on seawater intrusion should be more pronounced than that assuming homogeneity.

For water circulations, we found that locations and rates of seawater exchange are highly uncertain in heterogeneous coastal aquifers. This is significantly different from the results for homogeneous aquifers (Pu et al., 2020). We don’t recommend using equivalent homogeneous models to investigate water circulations in heterogeneous aquifers. This is also demonstrated in recent studies (Michael et al., 2016; Kreyns et al., 2020; Yu and Michael, 2022). The temperature modulates density-driven circulations and influences total seawater influxes. For heterogeneous aquifers, the increased dispersion induced by spatial heterogeneity weakens density gradients (Solórzano-Rivas et al., 2019). However, warmer freshwater or seawater could affect hydraulic conductivity and increase seawater circulation rates. This enhances the uncertainty of total groundwater discharge caused by complex salinity differences and geological heterogeneity. Seawater circulation induced by temperature change might be an important contributor to groundwater discharge of coastal aquifers (Guimond et al., 2022).

For natural systems, freshwater temperatures are relatively stable, whereas seawater temperatures can vary significantly (Locarnini et al., 2013; Benz et al., 2017). With a higher seawater temperature, two opposing seawater circulation cells induced by double diffusion of salt and heat are evidenced by Pu et al. (2020). A variety of system properties (e.g., permeability, thermal diffusivity, heterogeneity), together with temperature contrasts, could affect the two opposing circulation cells. In this study, we demonstrated that the anticlockwise seawater circulation cell mainly induced by heat diffusion was sensitive to and negatively correlated with the log-permeability variance of heterogeneous aquifer. For natural systems with spatial heterogeneity, FW-SW temperature contrasts in most coastal aquifers might not be such significant as we used in this study. Hence, two opposing circulation cells induced by double diffusion of salt and heat would be less pronounced. However, for heterogeneous aquifers, freshwater discharge zones under uniform permeability might switch into the zones of seawater circulations. This phenomenon is enhanced by the bigger permeability contrast (e.g., the larger log-permeability variance in this study). The mixing of land-sourced and sea-sourced chemicals could be strengthened in highly heterogeneous aquifers. Furthermore, freshwater and seawater temperatures would further change biogeochemical reactions in coastal aquifers (Cogswell and Heiss, 2021), influencing marine ecosystems (e.g., causing coastal eutrophication, algal blooms (Santos et al., 2021)).

Conclusion

We have examined the effects of temperature on salinity distributions and water circulations in heterogeneous aquifers through Monte-Carlo realizations of permeability fields. Based on our results, the following conclusions could be drawn:

(1) The averaged salinity distributions in heterogeneous aquifers are significantly affected by thermal regime, particularly for warmer freshwater. The random permeability field could increase the thermal impact of warmer freshwater and thus lower the lateral hydraulic gradient further, in comparison with the results from homogeneous aquifers. This causes the freshwater-seawater interface to intrude landward significantly in heterogeneous aquifers.

(2) Water exchange across the aquifer-sea interface is also changed with thermal regime of heterogeneous aquifers. Randomly-distributed permeability increases salt dispersion and decreases density gradients, influencing water circulations of heterogeneous aquifers. With warmer seawater, two opposing circulation cells induced by double diffusion of salt and heat are weakened by random permeability fields.

(3) Salinity distributions and water circulations in heterogeneous coastal aquifers are affected by variances and correlation lengths of permeability fields. The increased permeability field variance could enhance clockwise seawater circulations and inhibit double diffusion of salt and heat. The larger correlation length could facilitate the thermal effects on salinity distributions.

It should be pointed out that the simulation results were based on fixed boundary temperatures. In natural systems, the thermal effects on water flow and salinity distributions are much more complex (e.g., time-variant temperature conditions) than that considered in this study. Dynamic forcing conditions, such as tides and varying freshwater input, also play important roles in adjusting groundwater hydrodynamic processes of coastal aquifers (Kuan et al., 2019; Yu et al., 2019). The impacts of thermal forcing on other types of heterogeneous aquifers (e.g., layered aquifers and aquifers with conduit and fractures) also deserve future studies. Nevertheless, this study has shed light on the combined effects of temperature and aquifer heterogeneity on flow and salinity distributions, which should be considered in managing coastal groundwater resources.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

XY: methodology, formal analysis, data curation, conceptualization, writing—original draft. PX: conceptualization, methodology, supervision, writing—review and editing, LP: conceptualization, methodology, writing—review and editing.

Funding

This study was supported by National Natural Science Foundation of China (U2040204), Natural Science Foundation of Jiangsu Province (BK20200020) and Fundamental Research Funds for the Central Universities (B220202059). XY acknowledges Project funded by China Postdoctoral Science Foundation (BX20220099 and 2022M711026).

Acknowledgments

The authors thank the reviewers for their valuable comments, which significantly helped us to improve the manuscript quality.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2022.1002587/full#supplementary-material

References

Abarca, E., Carrera, J., Sanchez-Vila, X., and Dentz, M. (2007). Anisotropic dispersive Henry problem. Adv. Water Resour. 30, 913–926. doi:10.1016/j.advwatres.2006.08.005

CrossRef Full Text | Google Scholar

Abdoulhalik, A., and Ahmed, A. A. (2018). Transience of seawater intrusion and retreat in response to incremental water-level variations. Hydrol. Process. 32, 2721–2733. doi:10.1002/hyp.13214

CrossRef Full Text | Google Scholar

Anderson, M. P. (2005). Heat as a ground water tracer. Ground Water 43, 951–968. doi:10.1111/j.1745-6584.2005.00052.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bear, J. (1979). Hydraulics of groundwater. New York: McGraw-Hill.

Google Scholar

Befus, K. M., Cardenas, M. B., Erler, D. V., Santos, I. R., and Eyre, B. D. (2013). Heat transport dynamics at a sandy intertidal zone. Water Resour. Res. 49, 3770–3786. doi:10.1002/wrcr.20325

CrossRef Full Text | Google Scholar

Bellin, A., and Rubin, Y. (1996). HYDRO_GEN: A spatially distributed random field generator for correlated properties. Stoch. Hydrol. Hydraul. 10, 253–278. doi:10.1007/BF01581869

CrossRef Full Text | Google Scholar

Benz, S. A., Bayer, P., and Blum, P. (2017). Global patterns of shallow groundwater temperatures. Environ. Res. Lett. 12, 034005. doi:10.1088/1748-9326/aa5fb0

CrossRef Full Text | Google Scholar

Blanco-Coronas, A. M., Duque, C., Calvache, M. L., and López-Chicano, M. (2021). Temperature distribution in coastal aquifers: Insights from groundwater modeling and field data. J. Hydrology 603, 126912. doi:10.1016/j.jhydrol.2021.126912

CrossRef Full Text | Google Scholar

Burnett, W. C., Bokuniewicz, H., Huettel, M., Moore, W. S., and Taniguchi, M. (2003). Groundwater and pore water inputs to the coastal zone. Biogeochemistry 66, 3–33. doi:10.1023/B:BIOG.0000006066.21240.53

CrossRef Full Text | Google Scholar

Chang, S. W., and Clement, T. P. (2012). Experimental and numerical investigation of saltwater intrusion dynamics in flux-controlled groundwater systems. Water Resour. Res. 48, 2012WR012134. doi:10.1029/2012WR012134

CrossRef Full Text | Google Scholar

Cogswell, C., and Heiss, J. W. (2021). Climate and seasonal temperature controls on biogeochemical transformations in unconfined coastal aquifers. JGR. Biogeosciences 126, 1–21. doi:10.1029/2021JG006605

CrossRef Full Text | Google Scholar

Cooper, H. H. (1959). A hypothesis concerning the dynamic balance of fresh water and salt water in a coastal aquifer. J. Geophys. Res. 64, 461–467. doi:10.1029/JZ064i004p00461

CrossRef Full Text | Google Scholar

Croucher, A. E., and O’Sullivan, M. J. (1995). The Henry problem for saltwater intrusion. Water Resour. Res. 31, 1809–1814. doi:10.1029/95WR00431

CrossRef Full Text | Google Scholar

Engineering ToolBox (2003a). Specific heat of some common substances. Retrieved from https://www.engineeringtoolbox.com/specific-heat-capacity-d_391.html (Accessed March 10, 2021).

Google Scholar

Engineering ToolBox (2003b). Thermal conductivity of selected materials and gases. Retrieved from https://www.engineeringtoolbox.com/thermal-conductivity-d_429.html (Accessed March 10, 2021).

Google Scholar

Freeze, R. A. (1975). A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. Water Resour. Res. 11, 725–741. doi:10.1029/WR011i005p00725

CrossRef Full Text | Google Scholar

Gelhar, L. W. (1986). Stochastic subsurface hydrology from theory to applications. Water Resour. Res. 22, 135S–145S. doi:10.1029/WR022i09Sp0135S

CrossRef Full Text | Google Scholar

Geng, X., Boufadel, M. C., Rajaram, H., Cui, F., Lee, K., and An, C. (2020a). Numerical study of solute transport in heterogeneous beach aquifers subjected to tides. Water Resour. Res. 56, 1–20. doi:10.1029/2019WR026430

CrossRef Full Text | Google Scholar

Geng, X., Heiss, J. W., Michael, H. A., Boufadel, M. C., and Lee, K. (2020b). Groundwater flow and moisture dynamics in the swash zone: Effects of heterogeneous hydraulic conductivity and capillarity. Water Resour. Res. 56, 28401. doi:10.1029/2020WR028401

CrossRef Full Text | Google Scholar

Geng, X., and Michael, H. A. (2020). Preferential flow enhances pumping‐induced saltwater intrusion in volcanic aquifers. Water Resour. Res. 56, 1–15. doi:10.1029/2019WR026390

CrossRef Full Text | Google Scholar

Glover, R. E. (1959). The pattern of fresh‐water flow in a coastal aquifer. J. Geophys. Res. 64, 457–459. doi:10.1029/JZ064i004p00457

CrossRef Full Text | Google Scholar

Goswami, R. R., and Clement, T. P. (2007). Laboratory-scale investigation of saltwater intrusion dynamics. Water Resour. Res. 43, 11. doi:10.1029/2006wr005151

CrossRef Full Text | Google Scholar

Guimond, J. A., Mohammed, A. A., Walvoord, M. A., Bense, V. F., and Kurylyk, B. L. (2022). Sea-level rise and warming mediate coastal groundwater discharge in the Arctic. Environ. Res. Lett. 17, 045027. doi:10.1088/1748-9326/ac6085

CrossRef Full Text | Google Scholar

Henry, H. R. (1964). Effect of dispersion on salt encroachment in coastal aquifers in Sea water in coastal aquifers. U.S. Geological Survey Water-Supply. Paper 1613, C70–C84.

Google Scholar

Hughes, J. D., and Sanford, W. E. (2005). SUTRA-MS: A version of sutra modified to simulate heat and multiple-solute transport: U.S. Geological survey open-file report 2004-1207, 141. Retrieved from: https://water.usgs.gov/nrp/gwsoftware/SutraMS/SutraMS.html (Accessed on April 20, 2021).

Google Scholar

Hughes, J. D., Vacher, H. L., and Sanford, W. E. (2007). Three-dimensional flow in the Florida platform: Theoretical analysis of Kohout convection at its type locality. Geol. 35, 663. doi:10.1130/G23374A.1

CrossRef Full Text | Google Scholar

Hunt, A. G., Skinner, T. E., Ewing, R. P., and Ghanbarian-Alavijeh, B. (2011). Dispersion of solutes in porous media. Eur. Phys. J. B 80, 411–432. doi:10.1140/epjb/e2011-10805-y

CrossRef Full Text | Google Scholar

Kerrou, J., and Renard, P. (2010). A numerical analysis of dimensionality and heterogeneity effects on advective dispersive seawater intrusion processes. Hydrogeol. J. 18, 55–72. doi:10.1007/s10040-009-0533-0

CrossRef Full Text | Google Scholar

Ketabchi, H., and Jahangir, M. S. (2021). Influence of aquifer heterogeneity on sea level rise-induced seawater intrusion: A probabilistic approach. J. Contam. Hydrol. 236, 103753. doi:10.1016/j.jconhyd.2020.103753

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohout, F. A. (1965). Section of geological sciences: A hypothesis concerning cyclic flow of salt water related to geothermal heating in the floridan aquifer. Trans. N. Y. Acad. Sci. 28, 249–271. doi:10.1111/j.2164-0947.1965.tb02879.x

CrossRef Full Text | Google Scholar

Kreyns, P., Geng, X., and Michael, H. A. (2020). The influence of connected heterogeneity on groundwater flow and salinity distributions in coastal volcanic aquifers. J. Hydrology 586, 124863. doi:10.1016/j.jhydrol.2020.124863

CrossRef Full Text | Google Scholar

Kuan, W. K., Xin, P., Jin, G., Robinson, C. E., Gibbes, B., and Li, L. (2019). Combined effect of tides and varying inland groundwater input on flow and salinity distribution in unconfined coastal aquifers. Water Resour. Res. 55, 8864–8880. doi:10.1029/2018WR024492

CrossRef Full Text | Google Scholar

Locarnini, R. A., Mishonov, A. V., Antonov, J. I., Boyer, T. P., Garcia, H. E., Baranova, O. K., et al. (2013)., Silver Spring, MD: World ocean atlas 2013, Vol. 1: Temperature. Retrieved from https://www.nodc.noaa.gov/cgi-bin/OC5/woa13/woa13.pl?parameter=thttps://www.nodc.noaa.gov/cgi-bin/OC5/woa13/woa13.pl?parameter=t (accessed on June 8, 2021).

Google Scholar

Lu, C., Chen, Y., Zhang, C., and Luo, J. (2013). Steady-state freshwater-seawater mixing zone in stratified coastal aquifers. J. Hydrology 505, 24–34. doi:10.1016/j.jhydrol.2013.09.017

CrossRef Full Text | Google Scholar

Mahmoodzadeh, D., and Karamouz, M. (2019). Seawater intrusion in heterogeneous coastal aquifers under flooding events. J. Hydrol. X. 568, 1118–1130. doi:10.1016/j.jhydrol.2018.11.012

CrossRef Full Text | Google Scholar

Micallef, A., Person, M., Berndt, C., Bertoni, C., Cohen, D., Dugan, B., et al. (2021). Offshore freshened groundwater in continental margins. Rev. Geophys. 59, 1–54. doi:10.1029/2020RG000706

CrossRef Full Text | Google Scholar

Michael, H. A., Scott, K. C., Koneshloo, M., Yu, X., Khan, M. R., and Li, K. (2016). Geologic influence on groundwater salinity drives large seawater circulation through the continental shelf. Geophys. Res. Lett. 43, 10782–10791. doi:10.1002/2016GL070863

CrossRef Full Text | Google Scholar

Moore, W. S., and Joye, S. B. (2021). Saltwater intrusion and submarine groundwater discharge: Acceleration of biogeochemical reactions in changing coastal aquifers. Front. Earth Sci. 9, 1–14. doi:10.3389/feart.2021.600710

CrossRef Full Text | Google Scholar

Nguyen, T. T. M., Yu, X., Pu, L., Xin, P., Zhang, C., Barry, D. A., et al. (2020). Effects of temperature on tidally influenced coastal unconfined aquifers. Water Resour. Res. 56, 1–17. doi:10.1029/2019WR026660

CrossRef Full Text | Google Scholar

Paldor, A., Aharonov, E., and Katz, O. (2020). Thermo-haline circulations in subsea confined aquifers produce saline, steady-state deep submarine groundwater discharge. J. Hydrology 580, 124276. doi:10.1016/j.jhydrol.2019.124276

CrossRef Full Text | Google Scholar

Pool, M., Post, V. E. A., and Simmons, C. T. (2015). Effects of tidal fluctuations and spatial heterogeneity on mixing and spreading in spatially heterogeneous coastal aquifers. Water Resour. Res. 51, 1570–1585. doi:10.1002/2014WR016068

CrossRef Full Text | Google Scholar

Pu, L., Xin, P., Nguyen, T. T. M., Yu, X., Li, L., and Barry, D. A. (2020). Thermal effects on flow and salinity distributions in coastal confined aquifers. Water Resour. Res. 56, 1–17. doi:10.1029/2020WR027582

CrossRef Full Text | Google Scholar

Pu, L., Xin, P., Yu, X., Li, L., and Barry, D. A. (2021). Temperature of artificial freshwater recharge significantly affects salinity distributions in coastal confined aquifers. Adv. Water Resour. 156, 104020. doi:10.1016/j.advwatres.2021.104020

CrossRef Full Text | Google Scholar

Santos, I. R., Chen, X., Lecher, A. L., Sawyer, A. H., Moosdorf, N., Rodellas, V., et al. (2021). Submarine groundwater discharge impacts on coastal nutrient biogeochemistry. Nat. Rev. Earth Environ. 2, 307–323. doi:10.1038/s43017-021-00152-0

CrossRef Full Text | Google Scholar

Sebben, M. L., Werner, A. D., and Graf, T. (2015). Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured Henry problem. Adv. Water Resour. 85, 93–108. doi:10.1016/j.advwatres.2015.09.013

CrossRef Full Text | Google Scholar

Simpson, M. J., and Clement, T. P. (2004). Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. Water Resour. Res. 40, 1–12. doi:10.1029/2003WR002199

CrossRef Full Text | Google Scholar

Smith, A. J. (2004). Mixed convection and density-dependent seawater circulation in coastal aquifers. Water Resour. Res. 40, 16. doi:10.1029/2003wr002977

CrossRef Full Text | Google Scholar

Solórzano-Rivas, S. C., Werner, A. D., and Irvine, D. J. (2019). Dispersion effects on the freshwater–seawater interface in subsea aquifers. Adv. Water Resour. 130, 184–197. doi:10.1016/j.advwatres.2019.05.022

CrossRef Full Text | Google Scholar

van Lopik, J. H., Hartog, N., Zaadnoordijk, W. J., Cirkel, D. G., and Raoof, A. (2015). Salinization in a stratified aquifer induced by heat transfer from well casings. Adv. Water Resour. 86, 32–45. doi:10.1016/j.advwatres.2015.09.025

CrossRef Full Text | Google Scholar

Voss, C. I., and Provost, A. M. (2008). A model for saturated-unsaturated, variable-density ground-water flow with solute or energy transport. Reston, United States: U. S. geological survey water-resources investigations report 02-4231, 270.

Google Scholar

Werner, A. D., Bakker, M., Post, V. E. A., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B., et al. (2013). Seawater intrusion processes, investigation and management: Recent advances and future challenges. Adv. Water Resour. 51, 3–26. doi:10.1016/j.advwatres.2012.03.004

CrossRef Full Text | Google Scholar

Wilson, A. M. (2005). Fresh and saline groundwater discharge to the ocean: A regional perspective. Water Resour. Res. 41, 1–11. doi:10.1029/2004WR003399

CrossRef Full Text | Google Scholar

Yu, X., and Michael, H. A. (2022). Impacts of the scale of representation of heterogeneity on simulated salinity and saltwater circulation in coastal aquifers. Water Resour. Res. 58, 1–18. doi:10.1029/2020WR029523

CrossRef Full Text | Google Scholar

Yu, X., and Michael, H. A. (2019). Mechanisms, configuration typology, and vulnerability of pumping-induced seawater intrusion in heterogeneous aquifers. Adv. Water Resour. 128, 117–128. doi:10.1016/j.advwatres.2019.04.013

CrossRef Full Text | Google Scholar

Yu, X., Xin, P., Wang, S. S. J., Shen, C., and Li, L. (2019). Effects of multi-constituent tides on a subterranean estuary. Adv. Water Resour. 124, 53–67. doi:10.1016/j.advwatres.2018.12.006

CrossRef Full Text | Google Scholar