Nat Hazards DOI 10.1007/s11069-017-2860-0 ORIGINAL PAPER

Impacts of inland boundary conditions on modeling seawater intrusion in coastal aquifers due to sea-level rise Dong-mei Sun1 • Si-xiang Niu1 • Yong-ge Zang1

Received: 3 August 2016 / Accepted: 5 April 2017 Ó Springer Science+Business Media Dordrecht 2017

Abstract Inland boundary conditions have significant impacts on seawater intrusion (SWI) due to sea-level rise (SLR). This study investigated the impacts of three inland boundary conditions (flux-controlled (FC), head-controlled (HC) and general-head boundary (GHB)) on modeling SWI due to SLR in coastal aquifers. First, we used the TOUGH2/EOS7 software program to solve the standard Henry’s problem and compared the results with those from the literatures. Numerical simulations were performed to investigate the impacts of the three inland boundary conditions on SWI induced by SLR in confined coastal aquifers (standard Henry’s problem) and unconfined coastal aquifers (Henry’s problem extended to unconfined domains). The simulation results show that HC systems in both confined and unconfined coastal aquifers are remarkably vulnerable to the effects of SLR owing to the resulting decreased head difference, while FC systems are relatively insensitive to SLR because the equivalent freshwater head at the inland boundary can increase by a similar magnitude to the SLR. The GHB boundary is more realistic than the other two boundary conditions: the recharge decreases because of SLR depending on the aquifer properties and the head difference between the reference head and the equivalent freshwater head at the seaside boundary. Keywords Seawater intrusion  Sea-level rise  Inland boundary conditions  General-head boundary  Numerical simulation

1 Introduction Seawater intrusion (SWI) occurs in many coastal regions around the world and is affected by sea-level rise (SLR) due to global warming. Recently, there has been increasing interest in exploring the problem of SWI due to SLR, and inland boundary conditions have been

& Dong-mei Sun [email protected] 1

State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China

123

Nat Hazards

found to have a significant impact on SWI due to SLR (e.g., Werner and Simmons 2009; Carretero et al. 2013; Morgan et al. 2015). The two boundary conditions, the flux-controlled (FC) systems (e.g., Feseker 2007; Watson et al. 2010; Chang et al. 2011; Chang and Prabhakar 2012; Michael et al. 2013; Luoma and Okkonen 2014; Chesnaux 2015) and the head-controlled (HC) systems (e.g., Sherif and Singh 1999; Giambastiani et al. 2007; Lu and Werner 2013; Green and MacQuarrie 2014) have been extensively explored in previous studies. Specifically, Feseker (2007) conducted a two-dimensional density-dependent solute transport model to investigate the impact of SLR rate of 0.5 m per century on SWI under FC systems. Chang et al. (2011) used SEAWAT to evaluate the impact of SLR on SWI under FC systems in both confined and unconfined aquifers. Chesnaux (2015) used closed-form analytical solutions under FC systems to analyze the consequences of SLR on SWI in sloping coastal aquifers. Sherif and Singh (1999) observed that SWI in HC systems enhanced 9.0 km in the Nile Delta Aquifer response to a 500-mm rise in the Mediterranean, however, only increased 0.4 km with the same rise in sea level in the Bay of Bengal due to the different hydrogeologic conditions. Lu and Werner (2013) numerically investigated the impacts of SLR 1 m on the timescales of the intrusion and retreat of seawater in HC coastal aquifers. Green and MacQuarrie (2014) used a three-dimensional numerical model to evaluate the relative importance of the effects of climate change and groundwater extraction on SWI in HC coastal aquifers due to SLR 0.93 and 1.86 m in Atlantic Canada. Furthermore, some researchers compared the effects of these two inland boundary conditions on the SWI induced by SLR (e.g., Werner and Simmons 2009; Carneiro et al. 2010; Werner et al. 2012; Ataie-Ashtiani et al. 2013; Carretero et al. 2013; Morgan et al. 2015). Werner and Simmons (2009) adopted a simple analytical study to demonstrate that the unconfined coastal aquifer with a HC systems was remarkably more vulnerable to the impact of SLR, in comparison to that with a FC systems under otherwise the same conditions. The similar results were also obtained by Carretero et al. (2013). Werner et al. (2012) found that HC systems were more sensitive to SLR than FC systems in both confined and unconfined aquifers, based on an existing analytical solution for the steady-state position of a sharp interface. An exhaustive review on this topic was provided by Ketabchi et al. (2016). Although FC systems and HC systems have been widely used in previous studies, both have inherent limitations in practice because actual inland boundaries with SLR are neither FC nor HC systems but rather are indirect processes. To solve this problem, a general-head boundary (GHB), which is the simplest head-dependent flux boundary, has been introduced to represent a more realistic boundary. Dausman and Langevin (2005) completed SEAWAT simulations for the GHB coastal aquifer in Broward Country, Florida, and indicated that if the SLR becomes greater than 48 cm over the next century then several local well field would be contaminated. Oude Essink et al. (2010) showed numerically the impact of future SLR on SWI under GHB conditions in the Netherlands, and found that the SLR impact would be limited to areas within 10 km in coastal aquifers on account of the increased head of groundwater. Also, Lu et al. (2015) provided an analytical solution for SWI using this boundary condition. However, to our knowledge, no generalized investigation based on hypothetical and idealized conceptual models has been reported. Furthermore, the three inland boundary conditions have not been compared based on hypothetical numerical simulations to predict the effects of SLR on SWI. The main objective of this study is to investigate the impacts of three inland boundary conditions on SWI induced by SLR in both confined and unconfined coastal aquifers based on Henry’s problem (Henry 1964) and to compare the results with analytical solutions. The numerical simulations were conducted using the TOUGH2/EOS7 software program. EOS7

123

Nat Hazards

is a module in the TOUGH2 simulator that is used to simulate the two-phase flow of saline water and air (Pruess 1991; Pruess et al. 1999).

2 Inland boundary conditions The general conceptual models for confined and unconfined aquifers in coastal areas are shown in Fig. 1. The distance from the seaside boundary to the inland boundary is L [L], the thickness of the confined aquifer is B [L] and the height of the water table above the bottom of the aquifer is hf [L]. Hf [L] and qf [L2/T] are the actual head and actual freshwater recharge rates per unit width at the inland boundary of the domain, respectively, Hs is the height of the mean sea level above the bottom of the aquifer and Xt [L] is the distance of the intrusion from the seaside boundary. The homogeneous and isotropic aquifers are bounded on the bottom by horizontal impermeable strata regardless of the surface recharge and by a stationary seawater body at the seaside boundary. Three different inland boundary conditions are used: FC systems, in which the freshwater recharge qf is constant; HC systems, which maintain the actual equivalent freshwater head Hf; and a head-dependent flux boundary condition (GHB). GHB conditions are applied in a groundwater model to represent heads that are influenced by a large surface water body (e.g., wetlands, streams/rivers, reservoirs) with a known water elevation (i.e., Href ) outside the model domain. The flow that enters or leaves the model domain is proportional to the difference between the head at the GHB and the reference head Href , and it can be expressed by a linear equation (Harbaugh et al. 2000): q ¼ CðH  Href Þ

ð1Þ

where q [L2/T] is the inland freshwater recharge per unit width, C [L/T] is the hydraulic conductance and H [L] is the head at the GHB. Lu et al. (2015) derived the relationships between Hf and qf for both confined and unconfined conditions, which are written as follows: Confined:   KB ð1 þ aÞ B Hf  Hs þ ð2Þ qf ¼ L a 2a Unconfined: qf ¼

  K ð1 þ aÞ 2 Hf2  Hs 2L a

ð3Þ

in which K [L/T] is the hydraulic conductivity of the aquifer, and a is the density ratio, which is expressed as qf/(qs - qf) (where qf [ML-3] and qs [ML-3] are the freshwater density and seawater density, respectively) and is equal to 40. For GHB conditions in coastal systems, Hf and qf can be determined by combining Eqs. (1) and (2) for confined aquifers and Eqs. (1) and (3) for unconfined aquifers and letting H = Hf and q = qf: Confined: CL

Hf ¼ KB

Þ B Href þ ð1þa a Hs  2a CL KB þ 1

ð4Þ

123

Nat Hazards

Fig. 1 Conceptual models of coastal aquifers: a confined, b unconfined

qf ¼

Þ BC CHref  ð1þa a CHs þ 2a CL KB þ 1

ð5Þ

Unconfined: Hf ¼

123

CL þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Þ 2 2 C2 L2 þ ð1þa a K Hs þ 2CKLHref K

ð6Þ

Nat Hazards

0 qf ¼ C @Href 

CL þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 Þ 2 2 C2 L2 þ ð1þa a K Hs þ 2CKLHref A K

ð7Þ

3 Numerical simulations of a confined aquifer 3.1 Henry’s problem 3.1.1 Description of Henry’s problem In this study, a simplified and hypothetical confined coastal aquifer was used to investigate SWI into a confined coastal aquifer, which is known as Henry’s problem (Henry 1964). This problem was studied in a vertical slice of an aquifer that is perpendicular to the coast and parallel to the steady-state flow path. Figure 2 shows the vertical section of the simulated aquifer and the boundary conditions for Henry’s problem. The aquifer was assumed to be homogeneous and isotropic and was bounded above and below by impermeable strata. A constant freshwater flux was applied at the inland boundary (left boundary), while the seaside boundary (right boundary) was bounded by a stationary seawater body. Henry’s problem has been studied by many researchers (e.g., Pinder and Cooper 1970; Lee and Cheng 1974; Segol et al. 1975; Frind 1982; Huyakorn et al. 1987; Voss and Souza 1987; Oldenburg and Karsten 1995; Cheng et al. 1998; Liu et al. 2001; Simpson and Clement 2003; Rastogi et al. 2004; Karasaki et al. 2006; Abarca et al. 2007; Hamidi and Sabbagh-Yazdi 2008; Abd-Elhamid and Javadi 2011). The horizontal length of the model domain is L = 200 m, and the depth is B = 100 m. The inland freshwater recharge per unit width is Q = 6.6 9 10-1 m2/d. The saturated hydraulic conductivity is K = 1.0 m/d, the freshwater density is qw = 1000 kg/m3 and the

Fig. 2 Vertical section of the simulated aquifer and boundary conditions for Henry’s problem

123

Nat Hazards

seawater density is qs = 1025 kg/m3. The coefficient of dispersion is D = 6.6 9 10-2 m2/ d. Therefore, the three dimensionless parameters for the standard Henry’s problem are n ¼ BL ¼ 2:0, a ¼ KQ1 d ¼ 0:263 (where K1 ¼ K qsqqw ) and b ¼ D Q ¼ 0:1. w

3.1.2 Numerical simulation of the standard Henry’s problem The TOUGH2/EOS7 software program, which has been described in detail by Pruess (1991), was used to simulate the standard Henry’s problem. The model domain was discretized into 800 uniform hexahedral elements with Dx = Dz = 5.0 m (Fig. 2) and a unit width, i.e., Dy = 1.0 m. Table 1 summarizes the parameters that were used in the numerical simulation. It is worth noting that the constant diffusivity tortuosity model s0 sb ¼s0 Sb (where s0sb is the tortuosity of the medium, which is divided into a soildependent factor s0 and a saturation-dependent factor sb) is adopted; i.e., the effective diffusion coefficients in the saturated liquid phase are Djeff;l ¼ /dlj ¼6:6  102 m2/d. The TOUGH2/EOS7 software program uses four primary variables for the single liquidsaturated phase: the liquid phase pressure Pl, the brine mass fraction in the liquid phase Xb, the air mass fraction in the liquid phase Xal and the temperature T. The aquifer system was assumed to be isothermal (T = 25 °C). The freshwater mass flux per unit width was Fin = Q  qW = 7.639 9 10-3 kg/s, which was evenly assigned to all of the elements on the inland boundary as sink terms. Dirichlet boundary conditions (Pl ¼ Patm þ qs gð100  zÞ, Xb = 1.0 and Xal = 1.0 9 10-10, where Patm is the atmospheric pressure (1.013 9 105 Pa), and z[L] is the potential head) were applied at the seaside boundary. No flows were considered at the top and bottom boundaries. The initial conditions were Pl ¼ Patm þ qw gð100  zÞ, Xb = 1.0 and Xal = 1.0 9 10-10. In this situation, there was no capillary pressure, and the relative permeability was equal to 1.0. In this case, the inland boundary conditions formed an FC system. By running this transient problem for a total simulation time of 10,000 days, steadystate conditions were obtained after 5164 days. The 0.5 isochlor was chosen to represent the steady-state seawater wedge. Figure 3 shows the 0.5 isochlor distributions from this simulation and other previously published solutions by Frind (1982), Huyakorn et al. (1987), Cheng et al. (1998), Liu et al. (2001), Hamidi and Sabbagh-Yazdi (2008) and Abd-

Table 1 Parameters used in the numerical simulation

Parameter

Value

/: porosity

0.35

djl : coefficient of molecular diffusion in water and brine (m2/day)

18.8571 9 10-2

aL: longitudinal dispersivity (m)

0.0

aT: transverse dispersivity (m)

0.0 1.009 9 10-3

l: fluid viscosity (Pa s) 2

g: gravitational acceleration (m/s )

9.81

k: hydraulic conductivity (m/day)

1.0

K: intrinsic conductivity (m2)

1.18 9 10-12

3

6.6 9 10-1

Q: freshwater recharge (m /d) 3

123

qw: density of freshwater (kg/m )

1000

qs: density of seawater (kg/m3)

1025

Nat Hazards

Fig. 3 0.5 isochlor distributions under steady-state conditions

Elhamid and Javadi (2011). The results from this analysis are consistent with the previous solutions.

3.1.3 Numerical simulations of HC systems For HC systems, the equivalent head calculated by Eq. (2) is Hf = 102.57 m, which corresponds to qf = 0.66 m2/d and Hs = 100 m. Therefore, the inland boundary conditions of the HC systems were set as Pl = Patm þqw gð102:57  zÞ, Xb = 1.0 and Xal = 1.0 9 10-10. The steady-state conditions in the HC systems were reached after 7956 days. The simulation results for the HC systems are nearly identical to those for the FC systems. Figure 4 shows the 0.5 isochlor distributions for both the FC systems and HC systems under steady-state conditions, which indicate that the initial steady-state conditions for both systems are equivalent.

3.1.4 Numerical simulations of GHB conditions The values of qf = 0.66 m2/d and Hf = 102.57 m are fixed at the inland boundary for the FC and HC systems in the confined coastal aquifer, respectively. To produce the same steady-state conditions for GHB systems in the confined coastal aquifer, Href = 105 m was assumed. By letting Hf = 102.57 m and qf = 0.66 m2/d, C = 0.2716 m/s (Eq. (1)) and, therefore, the initial steady-state conditions prior to SLR for the GHB are the same as those for the FC and HC systems.

3.2 Effects of SLR in the confined coastal aquifer Section 3.1 showed that the initial steady-state conditions prior to SLR are equivalent with the three inland boundary conditions. To compare the impacts of SLR on SWI, different sea levels Hs of 100.25, 100.50, 100.75 and 101.0 m, which correspond to SLRs of 0.25, 0.50, 0.75 and 1.0 m, respectively, were applied to the seaside boundary. Despite the SLR, the freshwater flux qf = 0.66 m2/d and Hf = 102.57 m were fixed at the inland boundary in the FC and HC systems, respectively, while the value of qf varied with Hs in the GHB

123

Nat Hazards

Fig. 4 0.5 isochlor distributions for two systems under steady-state conditions

systems and was calculated with Eq. (5). The values of qf are 0.61, 0.57, 0.52 and 0.48 m2/ d for Hs values of 100.25, 100.50, 100.75 and 101.0 m, respectively. The inland boundary conditions for the GHB systems are similar to those of the FC systems except for the varying qf at the inland boundary. Figure 5 compares the steady-state 0.5 isochlor distributions for the three inland boundary conditions for SLRs of 0.25, 0.50, 0.75 and 1.0 m in the confined coastal aquifer, and Fig. 6 shows the 0.5 isochlor distributions after the SLR of 0.5 m for different inland boundary conditions. The SWI is relatively insensitive to SLR in the FC systems, but it is remarkably vulnerable to the SLR in the HC systems. The extent of the SWI in the GHB systems is between those of the FC and HC systems. To further explain these results, the equivalent freshwater head and freshwater recharge (per unit width) in the confined aquifer due to SLR with the three inland boundary conditions are compared in Table 2, in which IB [L] and SB [L] are the equivalent freshwater heads at the inland and seaside boundaries, respectively, and IB - SB [L] is the head difference between IB and SB. The equivalent freshwater heads at the inland freshwater flux boundary are obtained by averaging the heads of twenty elements located in the left column of the domain (the head in each element is calculated by h = z ? (Pl - Patm)/ (qwg)). The equivalent freshwater heads at the seaside boundary can be calculated with qffiffiffiffiffiffiffiffiffiffiffi h = Hs qs=q . Figure 7 compares the increases in water level at the inland boundary for w the three inland boundary conditions due to SLR in the confined aquifer. The magnitudes of the equivalent freshwater head difference and freshwater recharge (per unit width) for the three inland boundary conditions decrease with increasing sea level in the order FC [ GHB [ HC. Moreover, the differences between the three systems increase with increasing SLR. The head difference and recharge for FC systems remain nearly constant despite the SLR (in Table 2) because the equivalent freshwater head at the inland boundary also increases by a similar amount as the SLR at the seaside boundary (Fig. 7), while the GHB and HC systems experience gradual reductions, and the extent for the HC systems is greater than that for the GHB conditions.

123

Nat Hazards

Fig. 5 0.5 isochlor distributions in the confined coastal aquifer induced by different SLRs for different inland boundary conditions: a FC systems, b GHB and c HC systems

123

Nat Hazards

Fig. 6 0.5 isochlor distributions after an SLR of 0.5 m with the three inland boundary conditions Table 2 Comparison of the equivalent freshwater head and freshwater recharge (per unit width) in the confined coastal aquifer due to SLR for the three inland boundary conditions SLR 0.25 m

SLR 0.5 m

FC

GHB

HC

FC

GHB

HC

IB

102.80

102.72

102.57

103.06

102.90

102.57

SB

101.50

101.50

101.50

101.75

101.75

101.75

IB - SB

1.30

1.22

1.07

1.31

1.15

0.82

Recharge (10-3 kg/s)

7.64

7.12

5.88

7.64

6.60

4.24

SLR 0.75 m

SLR 1.0 m

FC

GHB

HC

FC

GHB

HC

IB

103.32

103.06

102.57

103.57

103.23

102.57

SB

102.00

102.00

102.00

102.25

102.25

102.25

IB - SB

1.32

1.06

0.57

1.32

0.98

0.32

Recharge (10-3 kg/s)

7.64

6.07

2.66

7.64

5.55

1.21

4 Numerical simulations of an unconfined coastal aquifer To explore the unconfined conditions, Henry’s problem is extended to unconfined domains that involve identical boundary conditions and parameters except for the top layer of the phreatic aquifer, which is increased by 5 m and has free surface that is discretized with the range of 0.25–1 m. As a result, the unconfined aquifer is 200 m long and has a total thickness of 105 m. To better examine the effects of SLR on the unconfined aquifer, the mesh was refined around the elevation of 100 m. Under unconfined conditions, the capillary pressure function is modeled using the van Genuchten function (van Genuchten 1980):

123

Nat Hazards

Fig. 7 Comparison of the increases in water level at the inland boundary for three inland boundary conditions due to SLR in the confined coastal aquifer

h i1k pc ¼ p0 ðS Þ1=k  1 ðpmax  pc  0Þ

ð8Þ

where p0 is the air entry pressure, k is a model parameter that relates to the degree of soil uniformity and S is the effective liquid phase saturation S* = (Sl - Slr)/(Sls - Slr), where Sl is the liquid saturation, Slr is the residual liquid saturation and Sls is the saturated liquid saturation. In this simulation, p0 = 1.373 9 103, k = 0.471, Sls = 0.921 and Slr = 0.146. The relative permeability function presented here is modeled using the van Genuchten– Mualem model (Mualem 1976; van Genuchten 1980): 8   k  2 ffi > < pffiffiffiffi  1=k  S 1  1  ðS Þ ðSl \Sls Þ krl ¼ ð9Þ > : 1 ðSl  Sls Þ The gas relative permeability krg can be chosen as one of the following two forms, the second of which is from Corey (1954): ( 1  krl ðSgr ¼ 0Þ krg ¼ ð10Þ 2 2 ^ ^ ðSgr [ 0Þ ð1  SÞ ð1  S Þ where S^ ¼ ðSl  Slr Þ=ð1  Slr  Sgr Þ, and Sgr is the residual gas saturation. In these simulations, Sgr = 0.079.

4.1 Numerical simulations to obtain the initial conditions The first step was to run the numerical model to obtain the steady-state initial conditions prior to SLR for the three inland boundary conditions. For the FC systems in the unconfined aquifer, the inland (0 B Z B 100 m) freshwater flux per unit width was Fin = 7.639 9 10-3 kg/s, which was the same as in the confined FC problem, and steadystate conditions were reached after 1217 days. To produce the same initial conditions with the HC systems in the unconfined aquifer, the constant head was Hf = 102.54 m (Eq. (3)), which corresponds to qf = 0.66 m2/d and Hs = 100 m, and the steady-state initial

123

Nat Hazards

Fig. 8 Steady-state 0.5 isochlor distributions for FC and HC systems in confined and unconfined aquifers

conditions were obtained after 7700 days. Similarly, for the initial conditions under GHB conditions, by letting qf = 0.66 m2/d, Hf = 102.54 m and Href = 105 m, the hydraulic conductance is C = 0.2681 m/s (Eq. (1)). Therefore, the initial conditions are equivalent for the three inland boundary conditions in the unconfined aquifer. Figure 8 compares the initial steady-state 0.5 isochlor distributions for the FC and HC systems in the confined and unconfined aquifers. The four curves are consistent with each other and nearly overlap, which indicates that the SWI behaviors under steady-state conditions are similar in confined and unconfined aquifers with both FC and HC systems. The initial steady-state conditions with the GHB are not shown because they are the same as in the FC systems.

4.2 Effects of SLR in the unconfined coastal aquifer In this section, the steady-state conditions for the three inland boundary conditions that were obtained in Sect. 4.1 were used as the initial conditions. To compare the effects of SLR on SWI in the unconfined aquifer, different sea levels Hs of 100.25, 100.50, 100.75 and 101.0 m, which correspond to SLRs of 0.25, 0.50, 0.75 and 1.0 m, respectively, were applied at the seaside boundary. Similarly, the freshwater flux qf = 0.66 m2/d and Hf = 102.54 m were fixed at the inland boundary despite the SLRs in the FC and HC systems, respectively, while the value of qf varied with Hs under the GHB conditions and were determined using Eq. (7). The values of qf are 0.62, 0.57, 0.53 and 0.48 m2/d for Hs values of 100.25, 100.50, 100.75 and 101.0 m, respectively. Figure 9 compares the steady-state 0.5 isochlor distributions in the unconfined aquifer with the three inland boundary conditions with SLRs of 0.25, 0.50, 0.75 and 1.0 m, and Fig. 10 shows the steady-state 0.5 isochlor distributions after a 0.5 m SLR with the different inland boundary conditions. The migrations of the SWI increase gradually under the three inland boundary conditions, which is similar to the results in the confined problems, and the migration with the GHB is between those of the other boundary conditions.

123

Nat Hazards

Fig. 9 0.5 isochlor distributions in the unconfined aquifer induced by different SLRs for the three inland boundary conditions: a FC systems, b GHB and c HC systems

123

Nat Hazards

Fig. 10 Steady-state 0.5 isochlor distributions after a 0.5 m SLR under different inland boundary conditions

The equivalent freshwater head and freshwater recharge (per unit width) in the unconfined aquifer due to SLR for the three inland boundary conditions are summarized in Table 3. Similar to the confined problems, the equivalent freshwater head difference and freshwater recharge (per unit width) decrease in the FC, GHB and HC systems with increasing sea level, and the effect increases with increasing SLR. Moreover, the decreases in the equivalent freshwater head difference and freshwater recharge (per unit width) in the HC systems are more pronounced than the other inland boundary conditions, the decrease with the GHB is between those of the other two boundary conditions, and there is little variation in the FC systems. Figure 11 compares the increases in water level at the inland boundary of the unconfined aquifer due to SLR for the three inland boundary conditions. The freshwater levels at the inland boundary in the FC systems are increased by a similar amount as the SLR and are more sensitive than under the GHB conditions, while they are fixed in the HC systems. Figure 12 shows the distributions of the phreatic level, which represents the free surface of the phreatic aquifer, for the three inland boundary conditions with different SLRs. The water levels at the inland boundary increased by approximately 0.9, 0.7 and 0 m from the initial state for an SLR of 1.0 m in the FC, GHB and HC systems, respectively.

5 Comparison with analytical solutions This section compares the results of the numerical simulations of GHB conditions in confined and unconfined coastal aquifers with analytical solutions. To compare the toe of the interface Xt using the numerical simulations with the analytical solution for the GHB conditions, a SLR of 1 m was simulated, and Href was varied between 103 and 120 m. The corresponding values of the interface toe Xt under the GHB conditions were calculated by Xt ¼ for the confined aquifer and

123

K 2 CB

þ LB 2aHref  2ð1 þ aÞHS þ B

ð11Þ

Nat Hazards Table 3 Comparison of the equivalent freshwater head and freshwater recharge (per unit width) in the unconfined coastal aquifer due to SLR for the three inland boundary conditions SLR 0.25 m

SLR 0.5 m

FC

GHB

HC

FC

GHB

HC

IB

102.76

102.62

102.54

102.99

102.81

102.54

SB

101.50

101.50

101.50

101.75

101.75

101.75

IB - SB

1.26

1.12

1.04

1.24

1.06

0.79

Recharge (10-3 kg/s)

7.64

7.13

6.78

7.64

6.62

5.23

SLR 0.75 m

SLR 1.0 m

FC

GHB

HC

FC

GHB

HC

IB

103.25

102.96

102.54

103.48

103.17

102.54

SB

102.00

102.00

102.00

102.25

102.25

102.25

IB - SB

1.25

0.96

0.54

1.23

0.92

0.29

Recharge (10-3 kg/s)

7.64

6.11

3.68

7.64

5.60

2.16

Fig. 11 Comparison of the increases in the water level at the inland boundary of the unconfined aquifer due to SLR for the three inland boundary conditions

Xt ¼

ð1 þ aÞK 2 Hs2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Þ 2 2 2a2 C Href K þ CL  C2 L2 þ ð1þa a K Hs þ 2CKLHref 

ð12Þ

for the unconfined aquifer (Lu et al. 2015). Figure 13 shows the analytical and numerical solutions for the distance of the interface toe (DXt) with the corresponding values of Href after a 1 m SLR in the confined and unconfined aquifers. There are significant differences between the analytical and numerical solutions in both the confined and unconfined aquifers. This is because of the different principles on which the solutions are based; specifically, the analytical solution assumes a sharp interface, while a transitional zone with a variable density model is assigned in the numerical simulation, so the numerical solution is represented by the steady-state 0.5

123

Nat Hazards

Fig. 12 Distributions of phreatic lines in the unconfined aquifer due to SLR for the three inland boundary conditions: a FC, b GHB and c HC

123

Nat Hazards

Fig. 13 Comparison of DXt with corresponding values of Href for a 1 m SLR

isochlor. In addition, the analytical solution can only provide the interface toe length Xt, but more information can be obtained from the numerical solution (e.g., the spatial and temporal variations of different isochlor distributions).

6 Discussion and conclusions This paper investigated the impacts of three inland boundary conditions on SWI induced by SLR in confined and unconfined coastal aquifers using the TOUGH2/EOS7 software program based on Henry’s problem. The simulation results indicate that the inland boundary conditions have significant impacts on SWI due to SLR. Both the confined and unconfined coastal aquifers are insensitive to the effects of SWI induced by SLR in FC systems because the equivalent freshwater head at the inland boundary can increase by a similar magnitude as the increase at the seaside boundary, which maintains the equivalent freshwater head difference between the inland and seaside boundaries. The aquifers are significantly more vulnerable to the SLR effect with the HC inland boundary owing to a reduction in freshwater recharge that is caused by the reduced equivalent freshwater head difference with the SLR. Moreover, the effect of SWI using the GHB is between those with the other inland boundary conditions and represents an intermediate boundary for investigating the effect of SLR on SWI. The GHB condition represents a more realistic boundary condition because the flux changes in response to the change in sea level. The key to the application of the GHB condition is to determine the appropriate value of the reference head. In general, the reference head can be determined from a large surface water body outside the domain with a known water elevation. If no large surface water body is present, the reference head may be identified by the phreatic level in the piedmont recharge area, in which the phreatic level is mainly affected by the local hydrogeological conditions, and the average value may remain nearly constant for many years. The numerical solutions were also compared with analytical solutions. There are significant differences between the solutions because of the different principles that are applied, and more information can be obtained from the numerical solutions. A better understanding of the impacts of different inland boundary conditions on SWI processes due

123

Nat Hazards

to SLR may help us to determine the proper inland boundary conditions in practical applications and to assess the potential effects of SLR on SWI. Acknowledgements This work was supported by the National Nature Science Foundations of China (Grant Nos. 51179118 and 51579170) and the Science Fund for Creative Research Groups of the National Natural Science Foundations of China (Grant No. 51321065), who provided financial support for the study design and advice in the data interpretation.

References Abarca E, Carrera J, Sa´nchez-Vila X, Dentz M (2007) Anisotropic dispersive Henry problem. Adv Water Resour 30(4):913–926 Abd-Elhamid HF, Javadi AA (2011) A density-dependant finite element model for analysis of saltwater intrusion in coastal aquifers. J Hydrol 401:259–271 Ataie-Ashtiani B, Werner AD, Simmons CT, Morgan LK, Lu C (2013) How important is the impact of landsurface inundation on seawater intrusion caused by sea-level rise? Hydrogeol J 21(7):1673–1677 Carneiro JF, Boughriba M, Correia A, Zarhloule Y, Rimi A, Houadi BE (2010) Evaluation of climate change effects in a coastal aquifer in Morocco using a density-dependent numerical model. Environ Earth Sci 61(2):241–252 Carretero S, Rapaglia J, Bokuniewicz H, Kruse E (2013) Impact of sea-level rise on saltwater intrusion length into the coastal aquifer, Partido de La Costa, Argentina. Cont Shelf Res 61–62:62–70 Chang SW, Prabhakar CT (2012) Experimental and numerical investigation of saltwater intrusion dynamics in flux-controlled groundwater systems. Water Resour Res 48(9):550–556 Chang SW, Clement TP, Simpson MJ, Lee KK (2011) Does sea-level rise have an impact on saltwater intrusion? Adv Water Resour 34(10):1283–1291 Cheng JR, Strobl RO, Yeh GT, Lin HC, Choi WH (1998) Modeling of 2D density-dependent flow and transport in the subsurface. Hydrol Eng 3(4):248–257 Chesnaux R (2015) Closed-form analytical solutions for assessing the consequences of sea-level rise on groundwater resources in sloping coastal aquifers. Hydrogeol J 23(7):1399–1413 Corey AT (1954) The interrelation between gas and oil relative permeability. Prod Mon 19(1):38–41 Dausman A, Langevin CD (2005) Movement of the saltwater interface in the surficial aquifer system in response to hydrologic stresses and water-management practices, Broward County, Florida. US Department of the Interior, US Geological Survey Feseker T (2007) Numerical studies on saltwater intrusion in a coastal aquifer in northwestern Germany. Hydrogeol J 15(2):267–279 Frind EO (1982) Simulation of long-term transient density-dependent transport in groundwater. Adv Water Resour 5:73–88 Giambastiani BMS, Antonellini M, Essink GHPO, Stuurman RJ (2007) Saltwater intrusion in the unconfined coastal aquifer of Ravenna (Italy): a numerical model. J Hydrol 340(1):91–104 Green NR, Macquarrie KTB (2014) An evaluation of the relative importance of the effects of climate change and groundwater extraction on seawater intrusion in coastal aquifers in Atlantic Canada. Hydrogeol J 22(3):609–623 Hamidi M, Sabbagh-Yazdi SR (2008) Modeling of 2D density-dependent flow and transport in porous media using finite volume method. Comput Fluids 37(8):1047–1055 Harbaugh AW, Banta ER, Hill MC, Mcdonald MG (2000) Modflow-2000, The U.S. Geological Survey Modular Ground-Water Model-User Guide to Modularization Concepts and the Ground-Water Flow Process. Open-file Report. U.S. Geological Survey 92:134 Henry HR (1964) Effects of dispersion on salt encroachment in coastal aquifers. Sea water in coastal aquifers 70 Huyakorn PS, Andersen PF, Mercer JW, White HO (1987) Saltwater intrusion in aquifers: development and testing of a three-dimensional finite element model. Water Resour Res 23:293–312 Karasaki K, Ito K, Maekawa K (2006) Simulation of salt water intrusion. Proceedings of the TOUGH Symposium 2006. Lawrence Berkeley National Laboratory, Berkeley, pp 15–17 Ketabchi H, Mahmoodzadeh D, Ataie-Ashtiani B, Simmons CT (2016) Sea-level rise impacts on seawater intrusion in coastal aquifers: review and integration. J Hydrol 535:235–255 Lee CH, Cheng TS (1974) On seawater encroachment in coastal aquifers. Water Resour Res 10(5):1039–1043

123

Nat Hazards Liu F, Turner I, Anh V (2001) A finite volume unstructured mesh method for modelling saltwater intrusion into aquifer systems. In: The first international conference and workshop on saltwater intrusion and coastal aquifers, monitoring, modeling, and management, Morocco Lu C, Werner AD (2013) Timescales of seawater intrusion and retreat. Adv Water Resour 59:39–51 Lu C, Xin P, Li L, Luo J (2015) Seawater intrusion in response to sea-level rise in a coastal aquifer with a general-head inland boundary. J Hydrol 522:135–140 Luoma S, Okkonen J (2014) Impacts of future climate change and baltic sea level rise on groundwater recharge, groundwater levels, and surface leakage in the Hanko Aquifer in Southern Finland. Water 6(12):3671–3700 Michael HA, Russoniello CJ, Byron LA (2013) Global assessment of vulnerability to sea-level rise in topography-limited and recharge-limited coastal groundwater systems. Water Resour Res 49(4):2228–2240 Morgan LK, Bakker M, Werner AD (2015) Occurrence of seawater intrusion overshoot. Water Resour Res 51(4):1989–1999 Mualem Y (1976) A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour Res 12(3):513–522 Oldenburg CM, Karsten P (1995) Dispersive transport dynamics in a strongly coupled groundwater-brine flow system. Water Resour Res 31:289–302 Oude Essink GHP, van Baaren ES, de Louw PGB (2010) Effects of climate change on coastal groundwater systems: a modeling study in the Netherlands. Water Resour Res 46:5613–5618 Pinder GF, Cooper HH (1970) A numerical technique for calculating the transient position of the saltwater front. Water Resour Res 6(3):875–882 Pruess K (1991) EOS7, an equation-of-state module for the TOUGH2 simulator for two-phase flow of saline water and air. Technical report LBL-31114, Lawrence Berkeley Laboratory, Berkeley, California Pruess K, Oldenburg CM, Moridis GJ (1999) TOUGH2 User’s Guide Version 2 Tough2 Users Guide Version. Lawrence Berkeley National Laboratory Rastogi AK, Choi GW, Ukarande SK (2004) Diffused interface model to prevent ingress of sea water in multi-layer coastal aquifers. Special Hydrol 42:1–31 Segol G, Pinder GF, Gray WG (1975) A Galerkin-finite element technique for calculating the transient position of the saltwater front. Water Resour Res 11(2):343–347 Sherif MM, Singh VP (1999) Effect of climate change on sea water intrusion in coastal aquifers. Hydrol Process 13(8):1277–1287 Simpson MJ, Clement TP (2003) Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models. Adv Water Resour 26:17–31 van Genuchten MTh (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc 44:892–898 Voss CI, Souza WR (1987) Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone. Water Resour Res 23(10):1851–1866 Watson TA, Werner AD, Simmons CT (2010) Transience of seawater intrusion in response to sea level rise. Water Resour Res 46:439–445 Werner AD, Simmons CT (2009) Impact of sea-level rise on sea water intrusion in coastal aquifers. GroundWater 47:197–204 Werner AD, Ward JD, Morgan LK, Simmons CT, Robinson NI, Teubner MD (2012) Vulnerability indicators of sea water intrusion. GroundWater 50(1):48–58

123