Evaluation of analytic solutions for steady interface flow where the aquifer extends below the sea

A computational approach is presented for steady Dupuit interface ﬂow where the aquifer extends below the sea. A detailed approach is outlined to determine the head at the coastline so that the solution below the leaky seabed may be combined with any type of steady Dupuit interface ﬂow in the aquifer below the land. The method allows for any inland boundary condition including speciﬁed head and speciﬁed ﬂux; cases of freshwater lenses caused by inﬁltration are also considered. The approach is implemented in a Python script and a Jupyter Notebook. Authors. creativecommons.org/licenses/by/4.0/).


Introduction
This Technical Note concerns analytic solutions for onedimensional steady Dupuit interface flow in coastal aquifers where the aquifer extends below the sea. The sea is separated from the aquifer by a leaky seabed. A variety of solutions have been published for steady interface flow where the aquifer extends below the sea (e.g., Edelman, 1972;Bruggeman, 1999;Kooi and Groen, 2001;Morgan et al., 2015). Sikkema and van Dam (1982) provided a detailed mathematical treatment, which was used by Bakker (2006) to derive a complete set of analytic solutions for the case where flow in the aquifer below the land is confined and uniform. Evaluation of the solution by Bakker (2006) is complicated. It requires determination of the type of flow (four types are distinguished), and when the tip of the interface reaches the end of the seabed, the solution requires evaluation of elliptic integrals and an iterative approach to determine parameters.
This Technical Note is, in part, a response to the recent calls for reproducibility in computational hydrology (Fienen and Bakker, 2016;Hutton et al., 2016;Barba, 2016), where a case is made that computational results cannot be reproduced or scrutinized when the source code is not available. Here, a cookbook recipe is provided for the evaluation of the part of the solution of Bakker (2006) in the aquifer below the sea. The solution below the sea can be coupled to any type of flow in the aquifer below the land, which may be simulated with, e.g., the Strack potential (Strack, 1976). The recipe is implemented in a Python computer program and combined with several options for the boundary conditions in the aquifer below the land. A Jupyter Notebook is developed to evaluate the position of the interface for a variety of cases. A Jupyter Notebook is an interactive document that integrates text, computer code, and results (Kluyver et al., 2016). The Python code and Jupyter Notebook are available from Bakker (2017).

Solution below the sea bottom
Consider one-dimensional steady Dupuit interface flow in a vertical cross-section (Fig. 1). The aquifer extends below the sea and the saltwater is at rest. The depth of the interface may be obtained from the head in the aquifer with the Ghijben-Herzberg equation.
Below the sea, the aquifer is bounded on top by a leaky layer separating the sea from the aquifer, so that flow is semiconfined. In cases where the leaky seabed is absent, the leaky layer represents the vertical resistance to flow of the aquifer (Anderson, 2005;Bakker, 2014). The leakage through the leaky layer is approximated as vertical and computed as where q z [L/T] is the vertical component of the specific discharge vector through the leaky layer, h is the freshwater head in the aquifer, h s is the freshwater head equivalent to the hydrostatic pressure in the saltwater at the top of the aquifer, and c [T] is the resistance to vertical flow of the leaky layer. The resistance c is computed from the thickness D and vertical hydraulic conductivity k v of the leaky seabed as c ¼ D=k v . In absence of a physical leaky layer, the resistance c represents the resistance to vertical flow of the aquifer (Bakker, 2014). The leaky layer may have a finite length L s or an infinite length. The hydraulic conductivity of the aquifer is k [L/T] and the thickness is H. The leakage factor k [L] is defined as The dimensionless density difference m s is defined as where q f and q s are the densities of freshwater and saltwater, respectively. The main parameters of the problem are summarized in Table 1. The flow in the aquifer below the land is not specified at this point. The discharge crossing the shoreline is called Q 0 [L 2 /T], but is often unknown prior to solving the problem. Separate solutions are used for flow in the aquifer below the sea and for flow in the aquifer below the land. First, solutions are presented for flow below the sea, which result in equations for the head in the aquifer at the shoreline in terms of Q 0 . A procedure to determine Q 0 from onshore boundary conditions is presented in a separate section. The shoreline is located at x ¼ 0 ( Fig. 1).
Equations are presented in terms of dimensionless variables. The dimensionless head / is defined as The dimensionless head as a function of the dimensionless coordinate x=k is governed by two dimensionless parameters, L s =k and l, where the latter is defined as Note that dimensionless parameter l is a combination of the discharge Q 0 crossing the shoreline, the aquifer parameters, and the dimensionless density difference m s .
Four different types of flow are distinguished depending on the position of the tip and the toe of the interface. The tip of the interface is the location where the interface touches the top of the aquifer, while the toe of the interface is the location where the interface touches the bottom of the aquifer (Fig. 1). For type I, the toe of the interface is in the aquifer below the land and the tip of the interface does not reach the end of the semi-confined layer (Fig. 2a). For type II, the toe of the interface is in the aquifer below the sea and the tip of the interface does not reach the end of the semiconfined layer (Fig. 2b). For type III, the toe of the interface is in the aquifer below the land, and the tip of the interface is at the end of the semi-confined layer (Fig. 2c). For type IV, the toe of the interface is in the aquifer below the sea and the tip of the interface is at the end of the semi-confined layer (Fig. 2d).
The type of flow is a function of L s =k and the dimensionless parameter l (Eq. (5)), which includes the discharge Q 0 crossing the shoreline. For example, when flow is of type I (Fig. 2a) and the discharge Q 0 increases, the toe of the interface moves towards the shoreline. If Q 0 is large enough, the toe will cross the shoreline (type II flow). The limiting case for which the toe is exactly at the shoreline is reached when l ¼ ffiffiffiffiffiffiffiffi 2=3 p , as derived by Bakker (2006). In the following, a cookbook recipe is presented to determine the type of flow. An outline of the cookbook recipe is given in Fig. 3. Equations are given for the dimensionless head / 0 at the shoreline and the length L of the outflow face for the different flow types. All equations are taken from Bakker (2006), where a detailed derivation is given. Following the recipe (Fig. 3), the first step is to compute the dimensionless parameter l. The flow is of type I if l < ffiffiffiffiffiffiffiffi 2=3 p (and the length of the semi-confining layer is long enough, which will be checked later), and the dimensionless head / 0 at the shoreline can be computed as The length of the outflow face L is The flow is of type II if l P ffiffiffiffiffiffiffiffi 2=3 p (and the length of the semiconfining layer is long enough), and / 0 and L can be computed as Next, it is checked whether the length of the semi-confining layer at the bottom of the sea is longer than the computed length of the outflow face L. If it is not longer, then the flow is of type III or IV. The calculations for type III and IV flow are more involved. First, the value of the parameter a t must be determined from the following equality: where fð/; aÞ ¼ ð3 À1=4 À 3 1=4 ÞFðh; jÞ þ 2 Á 3 1=4 Eðh; jÞ  where Fðh; jÞ and Eð/; jÞ are incomplete elliptic integrals of the first and second kind, respectively, and Eq. (11) is implicit in the parameter a t and needs to be solved using a numerical rootfinding routine.
The flow is of type III if l < ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1 þ a 3 t Þ=3 p , and / 0 may be computed as where a is the root of the equation ffiffiffiffiffiffiffiffiffiffiffi 3a=2 p ½fð/ 0 ; aÞ À fð0; aÞ þ L s =k ¼ 0 Eqs. (14) and (15) must be solved simultaneously using a numerical rootfinding routine.

Evaluation of elliptic integrals
Standard libraries can be used to evaluate the incomplete elliptic integrals in (12). Two notations are used for the incomplete elliptic integrals (e.g., DLMF, 2016, Section 19.1) and they are related as follows Fð/; jÞ ¼ Fð/jmÞ Eð/; jÞ ¼ Eð/jmÞ ð 19Þ The parameter j in the notation with the comma and the parameter m in the notation with the vertical bar are related by m ¼ j 2 . In the paper by Bakker (2006), the scipy package for Python was used to evaluate the elliptic integrals (Jones et al., 2001). In the scipy package, the latter form of the elliptic integrals, with the parameter m, is implemented. Unfortunately, Bakker (2006) did not realize the difference between the two notations, and evaluated the elliptic integrals by passing j rather than j 2 for the parameter m, resulting in incorrect figures for the type III and IV flow examples and the transition diagram (Figs. 5-8 in Bakker, 2006).

Implementation
In the previous, a cookbook recipe was presented to determine the dimensionless head / 0 at the shoreline in terms of the aquifer parameters and the dimensionless parameter l, which includes the discharge Q 0 crossing the shoreline. In case the flux at the shoreline is known (flux-specified condition), for example from a water balance, the solution for interface flow in the aquifer below the land can be obtained directly with, e.g., the Strack potential (Strack, 1976). The flow below the land can be confined flow or unconfined flow with or without an interface. The boundary conditions are that the head at the shoreline is equal to h 0 (obtained from / 0 ) and the other boundary conditions must be chosen such that the flux across the shoreline is equal to Q 0 . For many cases, the head is known at some point in the aquifer below the land (head-specified condition), but the flux towards the coast is unknown prior to solving the problem. For such cases, the solution is obtained in an iterative manner. A numerical   2. Uniform flow in a coastal aquifer that is confined below the land. The aquifer extends below the sea and is separated from the sea by a leaky layer (light grey). The interface is shown with a blue line for four types of flow: a) Type I flow with l ¼ 0:2, b) Type II flow with l ¼ 1:5, c) Type III flow with l ¼ 0:2 and Ls=k ¼ 0:8, and d) Type IV flow with l ¼ 1:5 and Ls=k ¼ 1:5. Vertical exaggeration is ten fold. rootfinding method is used to determine the value of Q 0 such that the solution meets the specified head at the specified point.
Both the flux-specified condition and the head-specified condition are implemented in a Jupyter Notebook. The presented recipe (Fig. 3) is implemented in Python so that the type of flow (type I, II, III, or IV) is determined automatically. The Jupyter Notebook can be used to simulate confined flow below the land (which may include an interface) where the discharge is uniform and equal to Q 0 . The Jupyter Notebook requires input of the hydraulic conductivity k, aquifer thickness H and absolute value of the head gradient G towards the coast upstream of the interface toe and computes Q 0 as Q 0 ¼ kHG.  Bakker (2006), but the results for type III and IV flow differ slightly from the results in Bakker (2006) where the elliptic integrals were evaluated incorrectly.
Finally, a solution is obtained for the same situation, but now the head is specified at x ¼ À1000 m. The shape of the interface is computed for three different specified heads, 0.25 m, 0.5 m, and 1 m (Fig. 4).

Unconfined flow with infiltration
As a final example, consider the case of unconfined interface flow below a long island where the flow is caused by a uniform infiltration rate N on the land (Fig. 5). The width of the island is W so that the outflow into the sea on either side of the island is Q 0 ¼ NW=2. The aquifer is so deep that the interface never reaches the bottom of the aquifer and the seabed is so long that the tip of the interface never reaches the end of the seabed. This means that flow is of Type 1. It also means that the flow does not depend on the thickness of the aquifer H. The head h 0 at the shoreline may be obtained from (6) by substituting (5) for l and canceling out terms, which gives where the elevation of the seabed is set to sea level (h s ¼ 0). Similarly, the length of the outflow face is obtained from (7) as The head and interface below the island may be obtained with the Strack potential. The discharge potential as a function of x (with the origin at the center of the island) may be written as where the relation between head and potential is The potential U 0 at the shoreline is obtained by substituting h 0 (20) for h in (23). As an example, the head and interface are computed for an island that is W ¼ 1 km wide. The hydraulic conductivity is k ¼ 10 m=d, the infiltration rate is N ¼ 0:001 m=d, the resistance of the leaky seabed is c ¼ 100 d, and the density of saltwater is q s ¼ 1025 kg/m 3 (Fig. 5).

Conclusions
Detailed guidelines were presented for the computation of steady interface flow below a leaky seabed. The solution below the seabed may be combined with a variety of solutions in the aquifer below the land. Examples were given for confined flow below the land with both a given uniform flux and a given inland head, and for unconfined interface flow below an elongated island with uniform areal infiltration. The presented computational approach is implemented in a Python script. A Jupyter Notebook is developed for interactive evaluation and visualization of the solutions for confined flow in the aquifer below the land including both flux-specified and head-specified conditions (Bakker, 2017). In addition. the Jupyter Notebook includes the presented example for an island aquifer system, which is expanded in the Notebook to include an impermeable base and the option to draw streamlines. Fig. 4. Uniform flow in a coastal aquifer that is confined below the land. The aquifer extends below the sea and is separated from the sea by a leaky layer (light grey). The interface position is shown for three different inland heads at x ¼ À1000 m. Aquifer thickness is 10 m. Vertical exaggeration is ten fold.