Effects of aquifer geometry on seawater intrusion in annulus segment island aquifers

Seawater intrusion in island aquifers was considered analytically, specifically for annulus segment aquifers (ASAs), i.e., aquifers that (in plan) have the shape of an annulus segment. Based on the Ghijben–Herzberg and hillslope-storage Boussinesq equations, analytical solutions were derived for steady-state seawater intrusion in ASAs, with a focus on the freshwater–seawater interface and its corresponding watertable elevation. Predictions of the analytical solutions compared well with experimental data, and so they were employed to investigate the effects of aquifer geometry on seawater intrusion in island aquifers. Three different ASA geometries were compared: convergent (smaller side is facing the lagoon, larger side is the internal no-flow boundary and flow converges towards the lagoon), rectangular and divergent (smaller side is the internal no-flow boundary, larger side is facing the sea and flow diverges towards the sea). Depending on the aquifer geometry, seawater intrusion was found to vary greatly, such that the assumption of a rectangular aquifer to model an ASA can lead to poor estimates of seawater intrusion. Other factors being equal, compared with rectangular aquifers, seawater intrusion is more extensive, and watertable elevation is lower in divergent aquifers, with the opposite tendency in convergent aquifers. Sensitivity analysis further indicated that the effects of aquifer geometry on seawater intrusion and watertable elevation vary with aquifer width and distance from the circle center to the inner arc (the lagoon boundary for convergent aquifers or the internal no-flow boundary for divergent aquifers). A larger aquifer width and distance from the circle center to the inner arc weaken the effects of aquifer geometry, and hence differences in predictions for the three geometries become less pronounced.

1 Introduction

Islands are extensively distributed throughout the world's oceans. Unfortunately, their groundwater resources are impacted by sea-level rise and increased demands. According to a recent estimate, there are approximately 65 million people living on oceanic islands where groundwater may be the only source of freshwater (Thomas et al., 2020). Fresh groundwater stored on oceanic islands is mainly from precipitation (usually in the form of a freshwater lens), and its availability varies due to different factors, e.g., island topography, rainfall patterns, tides, episodic storms and human activities (White and Falkland, 2010; Storlazzi et al., 2018). Seawater intrusion is thus an important issue due to its deleterious effect on oceanic island freshwater storage (e.g., Werner et al., 2017; Lu et al., 2019; Memari et al., 2020).

Figure 1 Island with an annulus segment in the Namu Atoll, Marshall Islands (© Google Earth 2021).

Over the past few decades, seawater intrusion in oceanic islands has been extensively investigated by field observations (e.g., Röper et al., 2013; Post et al., 2019), laboratory experiments (e.g., Stoeckl et al., 2015; Bedekar et al., 2019; Memari et al., 2020), numerical simulations (e.g., Lam, 1974; Gingerich et al., 2017; Liu and Tokunaga, 2019) and analytical solutions (e.g., Fetter, 1972; Ketabchi et al., 2014; Lu et al., 2019). Among these, analytical solutions are effective tools to assess the extent of seawater intrusion (i.e., the location of the freshwater–seawater interface), although they cannot incorporate complex factors (e.g., dispersive mixing and transient oceanic dynamics) (Werner et al., 2013). The advantages of analytical solutions are that they are computationally efficient, can be used as test cases for numerical models and can reveal the explicit relationships between parameters that influence seawater intrusion (e.g., Fetter, 1972; Ketabchi et al., 2014; Liu et al., 2014; Lu et al., 2019).

Based on the Dupuit–Forchheimer approximation (i.e., ignoring vertical flow) and the Ghijben–Herzberg equation (Drabbe and Badon Ghijben, 1889, English translation given by Post, 2018; Herzberg, 1901), Fetter (1972) presented analytical solutions describing the freshwater–seawater interface location and watertable elevation in a circular island. Bailey et al. (2010) further compared these single-layered analytical solutions with field measurements, indicating that the analytical solutions perform well in estimating the freshwater–seawater interface location and watertable elevation. Fetter's solutions formed the foundation for many subsequent analytical studies on seawater intrusion in island aquifers. Again, for a single layer, Chesnaux and Allen (2008) and Greskowiak et al. (2013) developed analytical solutions to predict the steady-state groundwater age distribution in freshwater lenses. In addition, using single-layered analytical solutions, Morgan and Werner (2014) proposed vulnerability indicators of freshwater lenses under sea-level rise and recharge change.

Since aquifers are usually heterogeneous, the single-layer analytical solutions were subsequently extended to two-layered island aquifers. Vacher (1988) derived solutions for the freshwater–seawater interface location and watertable elevation for infinite-strip islands composed of different layers. Dose et al. (2014) conducted laboratory experiments to validate and confirm the reliability of analytical solutions proposed by Fetter (1972) and Vacher (1988). Ketabchi et al. (2014) extended Fetter's analytical solutions to calculate the freshwater–seawater interface location and watertable elevation in two-layered circular islands subject to sea-level rise. Their results indicated that land-surface inundation caused by sea-level rise has a considerable impact on fresh groundwater lenses. Recently, Lu et al. (2019) derived analytical solutions for the freshwater–seawater interface location and watertable elevation for both strip and circular islands with two adjacent layers, i.e., a less permeable slice along the shoreline of an island and a more permeable zone inland.

All the abovementioned analytical solutions apply to either strip or circular islands. According to the classification of sand dunes developed by Stuyfzand (1993, 2017), there are different island layouts that should be considered, e.g., where the shape of the island is an annulus segment, instead of a strip or circular disk (Fig. 1). Annulus segment-shaped islands are found in various atolls (i.e., circular chains of islands surrounding a central lagoon) as found in the Pacific and Indian oceans (Werner et al., 2017; Duvat, 2019). Nevertheless, analytical solutions of seawater intrusion are not yet available for annulus segment aquifers (ASAs). In general, ASAs are conceptually treated as a 2D cross section, similar to strip islands (e.g., Ayers and Vacher, 1986; Underwood et al., 1992; Bailey et al., 2009; Werner et al., 2017). Evidently, topography plays an important role in groundwater flow and hence seawater intrusion (e.g., Zhang et al., 2016; Liu and Tokunaga, 2019). It remains unclear whether analytical solutions of seawater intrusion for strip islands are appropriate for ASAs. It is also unclear how island geometry affects the freshwater–seawater interface location and watertable elevation of ASAs.

In this study, analytical solutions are derived for steady-state seawater intrusion for ASAs, with a focus on the freshwater–seawater interface location and its corresponding watertable elevation. After comparing their predictions with experimental data (Memari et al., 2020), the analytical solutions are employed to investigate the effects of aquifer geometry on the freshwater–seawater interface location and watertable elevation in ASAs.

2 Conceptual model

Figure 2 shows the conceptual model of an ASA (a slice of an atoll island). The plan view of the model domain is represented as a sector (EFGH) with an angle θ (Fig. 2a). The sea (EF) and lagoon (HG) boundaries are located at L+L₀ (L) and L₀ (L) from the circle center, respectively. Since the longitudinal length is usually much longer than the lateral length for an atoll island (Werner et al., 2017), seawater intrusion from the lateral sides (EH and FG, Fig. 2a) is negligible in comparison to the longitudinal side, especially for the middle portion of an ASA. Therefore, EH and FG are treated as lateral no-flow boundaries. Note that treating the lateral sides as no-flow boundaries is often used in studies of freshwater lenses on atoll islands (e.g., Ayers and Vacher, 1986; Underwood et al., 1992; Bailey et al., 2009; Werner et al., 2017). The lateral vertical cross section of the model domain is conceptualized as a rectangle (ABCD) along the radial direction with dimensions of L (L) (width) × d (L) (height) (Fig. 2b and c). AD is the impermeable base, while BC is the land surface through which aquifer recharge flows.

Figure 2 Conceptual model of an annulus segment aquifer (a slice of an atoll island). (a) Plan view and (b, c) lateral vertical cross section with the saltwater interface tip (b) above the aquifer bed (single location) and (c) on the aquifer bed (two locations). In (a), the sea boundary is on EF, and the atoll lagoon boundary is on HG; in (b) and (c), AD is the impermeable base, and OO^* is the internal no-flow boundary.

Both the sea and lagoon water levels are set to H_s (L), which results in an internal no-flow boundary (water divide, where the slope of the watertable is zero) between the sea and lagoon (location of the z axis in Fig. 2b and c). The segment between the sea and the internal no-flow boundary is referred to as Unit 1, whereas the segment between the internal no-flow and lagoon boundaries is referred to as Unit 2 (Fig. 2). The widths of Units 1 and 2 are l₁ (L) and l₂ (L), respectively. In addition, the flow is asymmetrical in Units 1 and 2, with divergent flow (the aquifer length w (L) increases along the flow direction) in Unit 1 and convergent flow (w decreases along the flow direction) in Unit 2.

The r−z coordinate origin is placed at the intersection of the internal no-flow boundary and impermeable base, with the r axis pointing to the circle center (radial direction) and the z axis pointing vertically upward. Further, ϕ (L) is the watertable height, h (L) is the vertical distance between the watertable and the interface, h_s (L) is the vertical distance between the sea level and the interface and $h_{\mathrm{c}}=H_{\mathrm{s}}-h_{\mathrm{s}}$ (L) is the vertical distance from the impermeable base to the interface for a given r (Fig. 2b and c). Constant recharge into the saturated zone, N (L T⁻¹), is assumed. There are two possibilities for the interface tip (i.e., the location where the freshwater–seawater interface connects to the z axis or the bottom boundary): above the aquifer bed (Fig. 2b) or on the aquifer bed (Fig. 2c). The r coordinates of the interface tip in Units 1 and 2 are denoted as r_t1 (L) and r_t2 (L), respectively (Fig. 2c). Note that $r_{t\mathrm{1}}=r_{t\mathrm{2}}=\mathrm{0}$ when the interface tip is above the aquifer bed, as in Fig. 2b.

Consistent with previous studies (e.g., Ketabchi et al., 2014; Lu et al., 2016, 2019), the following assumptions are made: (1) steady-state flow, (2) sharp freshwater–seawater interface, (3) homogeneous and isotropic aquifer with a horizontal bottom, (4) aquifer recharge, due to rainfall infiltration, assumed to be stationary and uniform throughout the domain, (5) vertical flow in the saturated zone negligible (the Dupuit–Forchheimer approximation) and (6) the same velocity assumed on the arc (w) for a given radial distance r, leading to radial flow only. Based on this last assumption, the 3D flow problem can be simplified to 1D, making it possible to consider geometry effects analytically (Fan and Bras, 1998; Paniconi et al., 2003; Troch et al., 2003).

3 Analytical solutions

Under the abovementioned assumptions, groundwater flow in an ASA (Fig. 2) can be described as (Fan and Bras, 1998; Paniconi et al., 2003; Troch et al., 2003)

$$\begin{array}{}\text{(1)}& -{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(wq\right)+Nw=\mathrm{0},\end{array}$$

where q (L² T⁻¹) is the radial flux per unit length along the radial direction r (L). Equation (1) is a special case of the hillslope-storage Boussinesq equation proposed by Troch et al. (2003). Paniconi et al. (2003) validated the hillslope-storage Boussinesq equation by comparing it with a 3D Richards' equation model and found that predictions of the hillslope-storage Boussinesq equation matched those of the 3D model for seven different geometries well. For conciseness, readers are referred to Paniconi et al. (2003) for more details about the validation. Subsequently, the hillslope-storage Boussinesq equation was used for different analyses (Hilberts et al., 2005, 2007; Hazenberg et al., 2015, 2016; Kong et al., 2016; Luo et al., 2018), all of which focus on hillslope aquifers where the aquifer bottom is usually sloping. The hillslope-storage Boussinesq equation assumes that groundwater flow is parallel to the aquifer bottom (the Dupuit–Forchheimer approximation). Therefore, it can be applied to coastal unconfined aquifers where the aquifer bottom slope is usually mild (Lu et al., 2016).

According to Darcy's law and the Dupuit–Forchheimer approximation, the freshwater flux in the aquifer segment between the seaward boundary and interface tip can be calculated as (ϕ is independent of z)

$$\begin{array}{}\text{(2)}& q=-\underset{{\mathrm{h}}_{\mathrm{c}}}{\overset{\mathit{\varphi }}{\int }}{K}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}r}}\mathrm{d}z=-K_{\mathrm{s}}\left(\mathit{\varphi }-h_{\mathrm{c}}\right){\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}r}},\end{array}$$

where K_s (L T⁻¹) is the saturated hydraulic conductivity.

3.1 Interface tip above the aquifer bed

We first consider the situation where the interface tip is above the aquifer bed (Fig. 2b). In Unit 1 where $w=\mathit{\theta }(L_{\mathrm{0}}+l_{\mathrm{2}}-r)$, substituting Eq. (2) into Eq. (1) and then integrating it gives

$$\begin{array}{}\text{(3)}& \begin{array}{rl}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}& \left[{\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)}^{\mathrm{2}}-{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}\right]N\\ & =-\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)K_{\mathrm{s}}\left(\mathit{\varphi }-h_{\mathrm{c}}\right){\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}r}}.\end{array}\end{array}$$

According to the Ghijben–Herzberg equation, the vertical thickness of the freshwater zone (h) in the interface zone is given by

$$\begin{array}{}\text{(4)}& h=\mathit{\varphi }-h_{\mathrm{c}}=(\mathrm{1}+\mathit{\alpha })\left(\mathit{\varphi }-H_{\mathrm{s}}\right),\end{array}$$

where $\mathit{\alpha }={\mathit{\rho }}_{\mathrm{f}}/({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{f}})$ is the dimensionless density difference, and ρ_f [M L⁻³] and ρ_s [M L⁻³] are the freshwater and seawater densities, respectively. Substitution of Eq. (4) into Eq. (3) yields

$$\begin{array}{}\text{(5)}& \begin{array}{rl}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}& \left[{\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)}^{\mathrm{2}}-{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}\right]N\\ & =-K_{\mathrm{s}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)(\mathrm{1}+\mathit{\alpha })\left(\mathit{\varphi }-H_{\mathrm{s}}\right){\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}r}}.\end{array}\end{array}$$

Rearranging Eq. (5) produces

$$\begin{array}{}\text{(6)}& \begin{array}{rl}-{\displaystyle \frac{\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)N}{\mathrm{2}}}& +{\displaystyle \frac{N{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)}}\\ & =-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha })\left(\mathit{\varphi }-H_{\mathrm{s}}\right){\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}r}}.\end{array}\end{array}$$

Integrating Eq. (6) leads to

$$\begin{array}{}\text{(7)}& \begin{array}{rl}-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}& \ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)Nr\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N{r}^{\mathrm{2}}+C_{\mathrm{1}}=-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha }){\displaystyle \frac{{\left(\mathit{\varphi }-H_{\mathrm{s}}\right)}^{\mathrm{2}}}{\mathrm{2}}},\end{array}\end{array}$$

where C₁ is the integration constant that is determined by the sea boundary condition (i.e., $r=-l_{\mathrm{1}}$, ϕ=H_s),

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{C}_{\mathrm{1}}& ={\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}\ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}+l_{\mathrm{1}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)l_{\mathrm{1}}N\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{l}_{\mathrm{1}}^{\mathrm{2}}N.\end{array}\end{array}$$

The relation between h_s and ϕ is given by

$$\begin{array}{}\text{(9)}& h_{\mathrm{s}}=\mathit{\alpha }\left(\mathit{\varphi }-H_{\mathrm{s}}\right).\end{array}$$

Combining Eq. (7) with Eq. (9) and eliminating ϕ yields

$$\begin{array}{}\text{(10)}& \begin{array}{rl}-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}& \ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)Nr\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N{r}^{\mathrm{2}}+C_{\mathrm{1}}=-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha }){\displaystyle \frac{h_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}}}.\end{array}\end{array}$$

Equation (10) gives the freshwater–seawater interface location in Unit 1 once l₁ and l₂ are determined.

Equation (7) applies to Unit 2 by replacing C₁ with C₂,

$$\begin{array}{}\text{(11)}& \begin{array}{rl}-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}& \ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)Nr\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N{r}^{\mathrm{2}}+C_{\mathrm{2}}=-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha }){\displaystyle \frac{{\left(\mathit{\varphi }-H_{\mathrm{s}}\right)}^{\mathrm{2}}}{\mathrm{2}}},\end{array}\end{array}$$

where C₂ is chosen to satisfy the lagoon boundary condition (r=l₂, ϕ=H_s),

$$\begin{array}{}\text{(12)}& C_{\mathrm{2}}={\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}\ln \left(L_{\mathrm{0}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)l_{\mathrm{2}}N-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{l}_{\mathrm{2}}^{\mathrm{2}}N.\end{array}$$

Combining Eqs. (9) and (11) and eliminating ϕ leads to

$$\begin{array}{}\text{(13)}& \begin{array}{rl}-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}& \ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)Nr\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N{r}^{\mathrm{2}}+C_{\mathrm{2}}=-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha }){\displaystyle \frac{h_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}}}.\end{array}\end{array}$$

Equation (13) gives the freshwater–seawater interface location in Unit 2 once l₂ is determined. Since the sea level and lagoon water level are the same, an internal no-flow boundary exists between the sea and lagoon, i.e.,

$$\begin{array}{}\text{(14)}& r=\mathrm{0},{\left(h_{\mathrm{s}}\right)}_{\mathrm{unit}\mathrm{1}}={\left(h_{\mathrm{s}}\right)}_{\mathrm{unit}\mathrm{2}},\end{array}$$

where (h_s)_unit 1 and (h_s)_unit 2 represent h_s in Units 1 and 2, respectively.

Combining Eqs. (10), (13) and (14) leads to expressions for l₁ and l₂,

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{l}_{\mathrm{1}}=L+L_{\mathrm{0}}-\sqrt{{\displaystyle \frac{\mathrm{2}L{L}_{\mathrm{0}}+L^{\mathrm{2}}}{\mathrm{2}\ln \left(L+L_{\mathrm{0}}\right)-\mathrm{2}\ln \left(L_{\mathrm{0}}\right)}}}\text{(16)}& {\displaystyle }& {\displaystyle }{l}_{\mathrm{2}}=\sqrt{{\displaystyle \frac{\mathrm{2}L{L}_{\mathrm{0}}+L^{\mathrm{2}}}{\mathrm{2}\ln \left(L+L_{\mathrm{0}}\right)-\mathrm{2}\ln \left(L_{\mathrm{0}}\right)}}}-L_{\mathrm{0}}.\end{array}$$

As indicated by Eqs. (15) and (16), the internal no-flow boundary between the sea and lagoon only depends on L and L₀. For known l₁ and l₂, Eqs. (10) and (13) can be employed to predict the freshwater–seawater interface location in Units 1 and 2, respectively. Once the interface location is determined, h and ϕ are given by

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }h={\displaystyle \frac{\mathrm{1}+\mathit{\alpha }}{\mathit{\alpha }}}{h}_{\mathrm{s}}\text{(18)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }={\displaystyle \frac{h_{\mathrm{s}}}{\mathit{\alpha }}}+H_{\mathrm{s}}.\end{array}$$

3.2 Interface tip on the aquifer bed

When the interface tip is on the aquifer bed, the location of the internal no-flow boundary remains the same as for the interface tip above the aquifer bed. The freshwater–seawater interface for Units 1 and 2 can be determined by Eqs. (10) and (13), respectively. Then, from Eq. (17), h at the aquifer segment between the sea boundary and the interface tip is determined. To calculate h for the aquifer segment between the interface tip and the internal no-flow boundary, the r coordinate of the interface tip is found. At the interface tip of Unit 1 (r=r_t1),

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{s}}=H_{\mathrm{s}}\text{(20)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }={\displaystyle \frac{\mathrm{1}+\mathit{\alpha }}{\mathit{\alpha }}}{H}_{\mathrm{s}}.\end{array}$$

With Eqs. (10) and (20), r_t1 is given by

$$\begin{array}{}\text{(21)}& \begin{array}{rl}-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}& \ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r_{t\mathrm{1}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)N{r}_{t\mathrm{1}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N{r}_{t\mathrm{1}}^{\mathrm{2}}=-C_{\mathrm{1}}-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha }){\displaystyle \frac{H_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}}}.\end{array}\end{array}$$

Let

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }a={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N\text{(23)}& {\displaystyle }& {\displaystyle }b=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)N\text{(24)}& {\displaystyle }& {\displaystyle }c=-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}\end{array}$$

and

$$\begin{array}{}\text{(25)}& m=-C_{\mathrm{1}}-K_{\mathrm{s}}(\mathrm{1}+\mathit{\alpha }){\displaystyle \frac{H_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}}},\end{array}$$

then Eq. (21) becomes

$$\begin{array}{}\text{(26)}& a{r}_{t\mathrm{1}}^{\mathrm{2}}+b{r}_{t\mathrm{1}}+c\ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r_{t\mathrm{1}}\right)=m,\end{array}$$

which is solved by a root-finding method.

The freshwater discharge for the aquifer segment between the interface tip and the internal no-flow boundary is calculated as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}& \left[{\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)}^{\mathrm{2}}-{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}\right]N\\ & =-\left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)K_{\mathrm{s}}\mathit{\varphi }{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}r}}.\end{array}\end{array}$$

Repeating the steps from Eqs. (3) to (7) gives

$$\begin{array}{}\text{(28)}& \begin{array}{rl}-{\displaystyle \frac{{\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)}^{\mathrm{2}}N}{\mathrm{2}}}& \ln \left(L_{\mathrm{0}}+l_{\mathrm{2}}-r\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(L_{\mathrm{0}}+l_{\mathrm{2}}\right)Nr\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}N{r}^{\mathrm{2}}+C_{\mathrm{3}}=-{\displaystyle \frac{K_{\mathrm{s}}}{\mathrm{2}}}{\mathit{\varphi }}^{\mathrm{2}},\end{array}\end{array}$$

where C₃ is determined by substituting Eq. (20) into Eq. (28). Then, Eq. (28) can be adopted to calculate h for the segment between the interface tip and the internal no-flow boundary where h=ϕ.

Similarly, the r coordinate of the interface tip in Unit 2 (r_t2) is obtained by substituting Eq. (19) into Eq. (13). Then, the watertable (h) of the aquifer segment between the interface tip and the internal no-flow boundary for Unit 2 is computed by repeating the steps from Eqs. (21) to (28).

4 Results and discussion

4.1 Validation of the analytical solutions

The analytical solutions were validated by comparing their predictions with experimental data compiled from Memari et al. (2020), who reported experiments carried out using a 15^∘ radial tank. The tank contained three distinct chambers: internal no-flow boundary condition, porous medium and constant-head boundary condition (i.e., sea or lagoon). The internal no-flow and seaward boundaries were respectively located at 10 and 55.5 cm from the circle center, i.e., 45.5 cm from the internal no-flow boundary to the constant-head boundary along the radial direction. Note that the experimental tank corresponds to Unit 1 of the radial aquifer with l₁=45.5 cm and l₂=0, so the analytical results were calculated using Eqs. (10) and (26). The thicknesses of the porous medium and sea level were 28 and 25 cm, respectively, with $K_{\mathrm{s}}=\mathrm{1.23}\times {\mathrm{10}}^{-\mathrm{2}}$ m s⁻¹. The measured saltwater and freshwater densities were respectively 1.015 and 0.999 g mL⁻¹, leading to α=62. Two different recharge events with constant N, $\mathrm{2.46}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{1.08}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹, were considered in the experiments.

Figure 3 shows the comparison between analytical and experimental results of the freshwater–seawater interface for different recharge events. In general, the analytical solution predicts the freshwater–seawater interface well for both recharge events, despite there being some differences between the analytical results and the measurements, particularly in the zone near the constant-head boundary ($r=-\mathrm{45}$ cm). These deviations are likely due to assumptions made in the analytical solution, i.e., (i) a sharp freshwater–seawater interface, (ii) ignoring the effect of freshwater discharge and (iii) neglecting the vertical flow (the Dupuit–Forchheimer approximation).

Figure 3 Comparison between analytical and experimental (data compiled from Memari et al., 2020) results for the freshwater–seawater interface location for different recharge events. Note that the left and right sides are the sea and internal no-flow boundaries, respectively.

4.2 Effects of aquifer geometry on seawater intrusion

Previous studies showed that boundary conditions play a critical role in estimates of seawater intrusion (Werner and Simmons, 2009; Lu et al., 2016). Therefore, the internal no-flow boundary between the sea and lagoon was examined for various ASAs. As indicated by Eqs. (15) and (16), this internal no-flow boundary depends only on L and L₀. The values of l₁ and l₂ calculated respectively from Eqs. (15) and (16) are shown in Fig. 4 for three typical values of L (500, 1000 and 2000 m), with L₀ varying from 10² to 10⁶ m. In general, the internal no-flow boundary deviates from the middle of the ASA. When L₀ is less than 10⁵ m, l₁ is larger than l₂ for the three different values of L, indicating an internal no-flow boundary closer to the lagoon boundary. For example, taking L=2000 m and L₀=100 m leads to l₁=1240 m and l₂=760 m, with a deviation of 240 m (12 % of 2000 m) from the middle of the ASA. When L₀ exceeds 10⁵ m, however, the location of the internal no-flow boundary can be approximated as being at the middle of the ASA for all considered values of L. This is in contrast to strip and circular aquifers, where the internal no-flow boundary is always in the middle of the aquifer due to symmetry.

Figure 4 Widths of Unit 1 and Unit 2 versus L₀ for aquifers with different total width L.

Since the internal no-flow boundary location between the sea and lagoon deviates from the middle of the ASA, we expect aquifer geometry to play a significant role in controlling seawater intrusion. As mentioned previously, ASAs can be convergent (Unit 1) or divergent aquifers (Unit 2) where the extent of seawater intrusion may be different. However, for strip aquifers, both Units 1 and 2 are rectangular with the same extent of seawater intrusion. Therefore, three geometries were compared in this study: convergent, rectangular and divergent (Fig. 5). These geometries have been widely examined in hillslope hydrology regrading to the effects of aquifer geometry on runoff generation (Troch et al., 2003; Kong et al., 2016; Luo et al., 2018). To present the results more conveniently, we placed the r−z coordinate origin at the intersection of the constant-head boundary (sea or lagoon) and the impermeable base, with the r axis pointing horizontally to the internal no-flow boundary and the z axis vertically upward (Fig. 5). In addition, the distance between the constant-head boundary and the internal no-flow boundary (aquifer width) is denoted as L^* (Fig. 5), while the other parameters remain the same.

Figure 5 Three-dimensional view of (a) convergent (smaller side facing the lagoon), (b) rectangular and (c) divergent aquifers (larger side facing the sea) compared in this study. L^* represents the distance from the sea/lagoon to the internal no-flow boundary, i.e., l₁ or l₂ in Fig. 2. The internal no-flow boundary corresponds to the z axis in Fig. 2.

Following previous studies (e.g., Lu et al., 2016, 2019), different cases were selected to show the effects of aquifer geometry on seawater intrusion (Cases 1 and 2 in Table 1). According to Werner et al. (2017), the width of atoll islands generally varies from 100 to 1500 m along the radial direction. In order to focus on the effects of aquifer geometry on seawater intrusion, the same L^* and L₀ were assumed for the three aquifers, with L^* and L₀ equal to 1000 and 200 m, respectively. Note that L₀ is the distance from the circle center to the lagoon boundary for convergent aquifers, whereas it represents the distance from the circle center to internal no-flow boundary for divergent aquifers hereafter. The sand characteristics were the same as in the experiments of Memari et al. (2020). Two recharge events were considered (Cases 1 and 2, Table 1). The freshwater–seawater interface was calculated using the analytical solutions for the three different aquifers. Note that the Appendix presents analytical solutions for seawater intrusion in strip aquifers deduced from Lu et al. (2019).

Table 1 List of parameters use in different simulations.

^* The parameter is varied: the range of L₀ is from 200 to 6000 m, whereas the range of L^* is from 600 to 1600 m.

Figure 6 shows the freshwater–seawater interface calculated for Cases 1 and 2. As can be seen, the extent of seawater intrusion is noticeably different for the three aquifer geometries. For high recharge ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹), the interface tip is located at around 500 m for the divergent aquifer, which is about twice the value of the rectangular aquifer and 6 times the value for the convergent aquifer (Fig. 6a). When the recharge decreases to $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹, the interface tip moves further landward for the three aquifers as expected, but the difference between results is still great (Fig. 6b). The interface tip is displaced above the aquifer bed for both the rectangular and divergent aquifers, while it remains on the aquifer bed for the convergent aquifer. Regardless of the recharge rate, the most landward freshwater–seawater interface occurs in the divergent aquifer and vice versa for the convergent aquifer. This underlines that aquifer geometry plays a major role in controlling seawater intrusion, and hence it is necessary to account for aquifer geometry in analyses of seawater intrusion.

Figure 6 Freshwater–seawater interface predicted by analytical solutions for three different aquifers with (a) high and (b) low recharge (Cases 1 and 2 in Table 1). Note that r=1000 m is the internal no-flow boundary in Fig. 5.

Figure 7 Sensitivity of (a, b) the locations of the freshwater–seawater interface and (c, d) watertable to L₀ for convergent (a, c) and divergent (b, d) aquifers. The arrow in each plot shows the direction of increasing L₀ (values given in a, used to produce the different curves). Note that predictions for rectangular aquifers are independent of L₀.

4.3 Sensitivity analysis

A sensitivity analysis was conducted to investigate to what extent aquifer geometry affects seawater intrusion. Since we focus on the effects of aquifer geometry on the locations of the freshwater–seawater interface and watertable, values of L₀ and L^* were varied, with other parameters kept constant. When conducting the sensitivity analysis of L₀, L^* was fixed at 1000 m, which is a typical value for ASAs (Werner et al., 2017). Figure 7 shows the sensitivity of the locations of the freshwater–seawater interface and watertable to changes in L₀ (Case 3, Table 1). The freshwater–seawater interface and watertable elevation are independent of L₀ for rectangular aquifers (Appendix A). However, the freshwater–seawater interface and watertable elevation differ greatly when varying L₀ for both convergent and divergent aquifers, highlighting that L₀ plays an important role in affecting seawater intrusion. Specifically, as L₀ increases, the freshwater–seawater interface moves more landward (larger $r/L^{*}$, Fig. 7a), and its corresponding watertable elevation decreases (Fig. 7c) for convergent aquifers. In contrast, for divergent aquifers, increasing L₀ moves the freshwater–seawater interface more seaward (smaller $r/L^{*}$, Fig. 7b), and its corresponding watertable elevation increases (Fig. 7d). For a given L₀, divergent aquifers have the largest extent of seawater intrusion and the lowest watertable elevation, and the opposite is the case for convergent aquifers (Fig. 7).

Figure 8 Sensitivity of (a–c) the locations of the freshwater–seawater interface and (d–f) watertable to L^* for convergent (a, d), rectangular (b, e) and divergent (c, f) aquifers. The arrow in each plot points to the increase in L^* values used to construct each curve (values indicated in a).

Regardless of the freshwater–seawater interface and watertable elevation, the deviation between rectangular aquifers and divergent or convergent aquifers is significant when L₀ is less than 2000 m (Fig. 7). For example, the r coordinate of the interface tip (z=0) is 262 m for the rectangular aquifer at L₀=200 m, whereas it is 78 (31 % of that in the rectangular aquifer) and 500 m (191 % of that in the rectangular aquifer) for the convergent and divergent aquifers, respectively. As L₀ increases, the deviation between the three aquifers decreases. When L₀=2000 m, the r coordinate of the interface tip is 262, 209 (80 % of that in the rectangular aquifer) and 318 m (121 % of that in the rectangular aquifer) for the rectangular, convergent and divergent aquifers, respectively. As L₀ increases to 6000 m, the freshwater–seawater interface and watertable elevation of both convergent and divergent aquifers tend to those of rectangular aquifers; i.e., geometry effects decrease with increasing L₀. These results highlight the critical role played by the shape of aquifers. As a result, ignoring the aquifer geometry may lead to an inappropriate management strategy for groundwater resources in atoll islands.

The sensitivity of the freshwater–seawater interface and watertable elevation to L^* was investigated by varying L^* from 600 to 1600 m while fixing L₀ to 200 m (Case 4, Table 1). As shown in Fig. 8, contrary to the results for varying L₀, in this case the freshwater–seawater interface and watertable elevation in all three topographies are related to L^*. Again, the extent of seawater intrusion is greatest in divergent aquifers and least in convergent aquifers for a given L^*. When L^* increases, the freshwater–seawater interface moves seaward, and the watertable elevation increases, regardless of aquifer geometry; i.e., the seawater intrusion decreases (Fig. 8a–c). This is because the total freshwater flux increases with increasing L^*, leading to a higher hydraulic gradient and hence less seawater intrusion (Fig. 8d–f). Moreover, an increase in L^* reduces the differences in the seawater intrusion distance among the three geometries; i.e., the effects of aquifer geometry on seawater intrusion are more significant at small L^*. However, even at the maximum L^* considered (1600 m), the deviation between three aquifers remains significant: the r coordinate of the interface tip is about 148 m for the rectangular aquifer, whereas it is about 32 (22 % of that in the rectangular aquifer) and 278 m (188 % of that in the rectangular aquifer) for the convergent and divergent aquifers, respectively. Both L₀ and L^* can greatly impact seawater intrusion estimates for divergent and convergent aquifers, highlighting the necessity to include geometry effects in analytical solutions of seawater intrusion.

5 Conclusions

Based on the Ghijben–Herzberg and hillslope-storage Boussinesq equations, we derived analytical solutions of steady-state seawater intrusion for ASAs, with a focus on the freshwater–seawater interface and its corresponding watertable elevation as affected by recharge. After comparing with experimental data of Memari et al. (2020), the analytical solutions were employed to examine the effects of aquifer geometry on seawater intrusion in island aquifers. Three different shapes of island aquifer were compared: convergent, rectangular and divergent. The results lead to the following conclusions.

• The presented analytical solutions perform well in predicting the experimental freshwater–seawater interface, suggesting that these analytical solutions can predict seawater intrusion reasonably in different aquifer geometries.

• Island geometry plays a significant role in affecting the freshwater–seawater interface and watertable elevation. Other factors being equal, the extent of seawater intrusion is greatest in divergent aquifers and conversely least in convergent aquifers. In contrast, the watertable elevation is lowest in divergent aquifers and highest in convergent aquifers.

• The effects of aquifer geometry on seawater intrusion are dependent on the aquifer width and distance from the circle center to the internal no-flow boundary (Figs. 7 and 8). A larger aquifer width and distance from the circle center to the inner arc (the lagoon boundary for convergent aquifers or the internal no-flow boundary for divergent aquifers) weaken the role played by aquifer geometry and hence lead to a smaller deviation of the extent of seawater intrusion between the three topographies.

Real island aquifers are expected to exhibit more complexity than considered here; e.g., they will have more complex shapes and are subjected to transient flow conditions caused by tides, waves and groundwater pumping (Mantoglou, 2003; Pool and Carrera, 2011; Werner et al., 2013). In addition, since the experimental scale of Memari et al. (2020) is necessarily small, future experiments and field data are needed to further validate and facilitate the analytical solutions. Despite this, the new analytical solutions, validated against experiments, can be used as a tool for rapid estimation of seawater intrusion in ASAs once known island geometry and corresponding soil properties are given.

Appendix A: Analytical solutions for rectangular aquifers

For rectangular aquifers, the seawater intrusion in Unit 1 is identical to that in Unit 2 because of symmetry. With the interface tip on the aquifer bed, analytical solutions for the freshwater–seawater interface (h_s), watertable elevation (h), and r coordinate of the interface tip in Unit 2 (r_t2) can be respectively written as (Lu et al., 2019)

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{s}}=\mathit{\alpha }\sqrt{{\displaystyle \frac{N}{(\mathrm{1}+\mathit{\alpha })K_{\mathrm{s}}}}\left({\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{4}}}-r^{\mathrm{2}}\right)}\text{(A2)}& {\displaystyle }& {\displaystyle }h=\left\{\begin{array}{ll}\sqrt{\frac{N}{K_{\mathrm{s}}}\left(r_{t\mathrm{2}}^{\mathrm{2}}-r^{\mathrm{2}}\right)+\left(\frac{H_{\mathrm{s}}}{\mathit{\alpha }}+H_{\mathrm{s}}\right)}& \mathrm{0}\le r\le r_{t\mathrm{2}}\\ \sqrt{\frac{N}{(\mathrm{1}+\mathit{\alpha })K_{\mathrm{s}}}\left(\frac{L^{\mathrm{2}}}{\mathrm{4}}-r^{\mathrm{2}}\right)}+H_{\mathrm{s}}& r_{t\mathrm{2}}<r\le \frac{L}{\mathrm{2}}\end{array}\right.\text{(A3)}& {\displaystyle }& {\displaystyle }{r}_{t\mathrm{2}}=\sqrt{{\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{4}}}-{\displaystyle \frac{(\mathrm{1}+\mathit{\alpha })K_{\mathrm{s}}}{N}}\left({\displaystyle \frac{H_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}}}\right)}.\end{array}$$

When the interface tip is above the aquifer bed, the analytical solution for the freshwater–seawater interface location and watertable elevation in Unit 2 are the same as Eqs. (A1) and (A2), respectively.