The Interface Between Fresh and Salt Groundwater in Horizontal Aquifers: The Dupuit–Forchheimer Approximation Revisited

We analyze the motion of a sharp interface between fresh and salt groundwater in horizontal, confined aquifers of infinite extend. The analysis is based on earlier results of De Josselin de Jong (Proc Euromech 143:75–82, 1981). Parameterizing the height of the interface along the horizontal base of the aquifer and assuming the validity of the Dupuit–Forchheimer approximation in both the fresh and saltwater, he derived an approximate interface motion equation. This equation is a nonlinear doubly degenerate diffusion equation in terms of the height of the interface. In that paper, he also developed a stream function-based formulation for the dynamics of a two-fluid interface. By replacing the two fluids by one hypothetical fluid, with a distribution of vortices along the interface, the exact discharge field throughout the flow domain can be determined. Starting point for our analysis is the stream function formulation. We derive an exact integro-differential equation for the movement of the interface. We show that the pointwise differential terms are identical to the approximate Dupuit–Forchheimer interface motion equation as derived by De Josselin de Jong. We analyze (mathematical) properties of the additional integral term in the exact interface motion formulation to validate the approximate Dupuit–Forchheimer interface motion equation. We also consider the case of flat interfaces, and we study the behavior of the toe of the interface. In particular, we give a criterion for finite or infinite speed of propagation.


Introduction
The study of simultaneous movement of fresh and salt groundwater is generally motivated by practical environmental problems such as seawater intrusion in coastal aquifers, upconing of saltwater under freshwater supply wells and seepage of brackish water in polder areas. An important aspect in such problems is the prediction of the position, and movement of the diffusive/dispersive mixing zone between fresh and salt groundwater in aquifers is required. The interest in fresh-salt groundwater problems dates back to the late nineteenth century, when Ghyben (1888) and Herzberg (1901) independently estimated the thickness of a fresh water lens in a coastal aquifer by considering a stationary fresh/salt interface, see also Verruijt (1980).
Field observations show that in many cases of practical interest the width of the mixing zone between the fluids is small compared to the vertical and horizontal extensions of the aquifer. This allows to disregard the mixing zone, resulting in a sharp interface model. Though relatively straightforward to formulate, fresh-salt interface problems pose many serious mathematical difficulties. The time-dependent case for incompressible fluids yields two-phase free boundary problems for which only partial existence and uniqueness results are known, see, e.g., Robert (1985, 1986) and Chan Hong and Hilhorst (1987). One of the main difficulties concerns the regularity of the free boundary. A detailed numerical approach was given by Chan Hong et al. (1989). Many saltwater intrusion problems can be considered stationary. Examples are given in Bear and Verruijt (1987) and by Dagan (2007, 2008). For stationary problems, rigorous mathematical results are known, see Alt (1990, 1993) and [17]. Choquet et al. (2015) considers an interesting approach for the time-dependent interface problem, see also Jazar and Monneau (2014). There is extensive engineering literature available. For this, we refer to, e.g., Dawson and Adnan (2015), Oude Essink (2001) and the books of Bear (1979), Bear et al. (1999) and Bear and Verruijt (1987). Recently, there is renewed interest in interface problems related to CO 2 -storage, see Nordbotten and Celia (2012).
Probably the most important simplification in groundwater flow is the so-called Dupuit-Forchheimer approximation, to which we will refer in this paper as the 'Dupuit approximation,' see Dupuit (1863) and Forchheimer (1886). The Dupuit approximation, originally formulated for groundwater in an unconfined aquifer with constant density, assumes that the pressure distribution below the groundwater table is hydrostatic, implying that the horizontal component of the specific discharge vector is constant over the depth. The Dupuit approximation reduces the complexity of a hydrogeological problem by reducing the number of space dimensions of the problem by one. The Dupuit formulation of a problem can also be the result of averaging of certain variables in the vertical direction, see, e.g., Strack et al. (2006), Huppert and Woods (1995), Bakker (1989), Strack and Bakker (1995), Danilov and Katz (1980) and in particular Yortsos (1995). More recently, sharp interface models in conjunction with the Dupuit approximation were studied in the framework of geological storage of CO 2 , see, e.g., Gasda et al. (2011), Nordbotten andDahle (2011), and for a rather complete overview the book by Nordbotten and Celia (2012).
Inspired by the unpublished work of the Dutch engineer J.H. Edelman, De Josselin de Jong (1960Jong ( , 1977Jong ( , 1969) developed a mathematical formulation for the dynamics of a two-fluid interface, using the stream function as flow variable. In this formulation, one replaces the two different homogeneous fluids by one hypothetical fluid with a distribution of vortices along the interface. The strength of the vortices, resulting from horizontal density gradients, is chosen to obtain the correct discharge field throughout the flow domain. Note that Kidder (1956) addressed a similar problem in an aquifer of finite length using the velocity potential and the method of images to derive a series expansion solution for the evolution of the interface.
In his paper, De Josselin de Jong (1981) [or in the collected papers of De Josselin de Jong (2006)] applied the Dupuit approximation to a sharp fresh-salt interface in an aquifer of infinite horizontal extent. He assumed that the horizontal component of the specific discharge vector q x (x, z, t) is constant with respect to the vertical position z in each fluid, with having a jump at the interface. In other words, he assumed that where q x f (x, t) and q xs (x, t) denote, respectively, the x-component of the specific discharge above (fresh) and below (salt) the interface. By parameterizing the interface along the horizontal coordinate and applying (1), he was able to derive an equation describing the movement of the fresh-salt interface. This equation is a nonlinear doubly degenerate diffusion equation, of which the mathematical properties have been studied in great detail, see, e.g., Van Duijn and Zhang (1988) and Zhang (1988). Note that we explicitly do not consider non-horizontal aquifers. For sake of completeness, we refer to the books of Charny (1963) and Charny et al. (1968) that cover dipping aquifers. Starting point for this paper is the exact stream function formulation for the flow field in the aquifer as presented in De Josselin de Jong (1981). Using this result, we are able to derive an exact integro-differential equation for the movement of the fresh-salt interface. It is shown that the pointwise (differential) part of the equation is identical to the interface motion equation as derived in De Josselin de Jong (1981) (and which was merely based on the physical intuition of De Josselin de Jong). The integral part of the exact equation can be regarded as the correction term to the approximate Dupuit equation. Our aim is to study the behavior of the integral term in order to give a criterion for the validity of the Dupuit approximation.
In Sect. 2, we formulate the flow problem and give an expression for the specific discharge in terms of an integral along the fresh-slat interface. Next in Sect. 3, we recall the expression for the Dupuit discharge and the related interface motion equation. In Sect. 4, we return to the interface integral and use integration by parts to recover the approximate Dupuit expression as well as a correction term that contains an integral involving a weight function and the second-order spatial derivative of the interface. Then, in Sect. 5, we present a criterion for the validity of the Dupuit approximation by estimating the correction term. The case of flat interfaces, including the explicit Dupuit similarity solution as a rotating interface, is considered in Sect. 6. Finally, Sect. 6.1 summarizes the main conclusion of this work.

Problem Formulation
Let the strip (2) denote a vertical cross section of a horizontally extended aquifer of constant thickness H , which is bounded above and below by impermeable confining layers. The vertical coordinate z is pointing upwards, i.e., opposite to the gravity vector. The aquifer consists of a homogeneous and isotropic porous medium, with intrinsic permeability κ and porosity φ , occupied by fresh and salt groundwater. We assume that the scale (H ) of the problem is sufficiently large, so that an abrupt change in specific weight γ , from freshwater with γ f to saltwater with γ s , 123 can be assumed. Typical values for seawater intrusion problems are: γ f = 1000 N/m 3 and γ s = 1025 N/m 3 . Throughout this paper, we consider the stable case where the heavier saltwater is always below the lighter freshwater. This implies that there exists an interface S between the fluids that can be parameterized as a function of the horizontal coordinate x. In other words, there exists a function ζ :R → [0, H ] such that For the specific weight γ , this implies or for any (x, z) ∈ Ω, where H denotes the Heaviside function, Finally, we assume that the fluids are incompressible and that they have the same viscosity μ, see also De Josselin de Jong (1960) and later Verruijt (1980). This generally holds in natural circumstances.

Flow Equations
Concerning the flow in Ω, we assume here without loss of generality that the fluids are at rest when |x| → ∞. Any superimposed flow can easily be incorporated. The simplifying assumptions lead to the following equations: Continuity div q = 0 in Ω; Momentum balance (Darcy's law) μ κ q + grad p + γ e z = 0 in Ω; Boundary conditions q · n = 0 on ∂Ω and |q| → 0 as x → ±∞ for z ∈ (0, H ).
Here, q denotes the specific discharge vector (with components q x and q z ), p the fluid pressure, e z the unit vector in the positive z-direction, n the outward unit vector normal to the boundary ∂Ω and | · | the Euclidian norm. Equation (7) is satisfied if a stream function Ψ :Ω → R exists such that Taking the two-dimensional curl of Darcy's law (8) and substituting (10) into the result gives the fundamental equation [see also the pioneering paper of De Josselin de Jong (1960)] Using expression (5) in this equation yields Here x ζ (x s ) = (x s , ζ(x s )) denotes a point at the interface S, x = (x, z) a point in the strip Ω and δ is the Dirac delta function in R 2 . Equation (12) tells us that the flow (stream function) induced by the fresh-salt interface is equivalent to the flow induced by vortices with strength Γ dζ dx (x) along the interface. This observation was one of the main results of De Josselin de Jong (1960). From the no-flow conditions (9), we deduce for Ψ the boundary conditions The unique, bounded solution of problem (12), (14) is given by the expression (see Appendix A for details and references) wherex = π x/H ,z = π z/H andx s = π x s /H . Hence, for a given interface ζ = ζ(x), the induced specific discharge components q x , q z are obtained by differentiating (15) according to (10). Expressing the cosine and hyperbolic cosine in terms of exponentials and introducing complex notation one finds for any point (x, z) ∈ Ω, provided (x, z) / ∈ S, i.e., not chosen on the interface. In (16) and in other integrals to come, ζ always denotes ζ(x s ). To obtain the specific discharge components at a point of the interface, say at (x p , ζ(x p )), we first calculate their value at some point (x p , z p ), with z p > ζ(x p ) or z p < ζ(x p ), and then pass to the limit An expression similar to (16) was used by De Josselin de Jong (1981) to compute the discharge for a specific case: a piecewise linear interface. Here, expression (16) is the starting point of the analysis.

Dupuit Approximation
A first main result of this paper is that (16) can be written as the sum of local terms and a new integral expression containing the second-order derivative of ζ . We show this in Sect. 4 by applying integration by parts to (16). For the horizontal discharge q x , we find, see also (42), and where Int(x, z; ζ ) denotes the integral given by (44) and ζ x = dζ dx . Later, when introducing time dependence, we have ζ = ζ(x, t). Then ζ x denotes the spatial derivative with respect to x: ζ x = ∂ζ ∂ x . For q z (x, z) a similar expression is given in (43). The local term in (17) is the Dupuit approximation or Dupuit discharge. Note that it does not depend on the vertical coordinate z in each fluid. This was used by De Josselin de Jong (1981). In that paper, he assumed q x to be constant in z in each fluid, as formulated in (1). Using this assumption, he derived (18) from the following argument. The fluid balance The shear flow at the interface (as a result of pressure continuity in Darcy's law) implies Combining (20) and (21) gives the Dupuit expression (18) with the obvious notation q f x = q D x (z > ζ ) and q sx = q D x (z < ζ ).

Movement of the Interface
So far, we considered stationary discharges resulting from a stationary fresh-salt distribution: ζ = ζ(x) only. We now introduce the time dependence to determine the evolution of the interface, starting from a given initial (say at t = 0) situation. From now on, we write ζ = ζ(x, t) and consequently all derivatives become partial derivatives. Let denote the total discharge of saltwater in the positive x-direction. Inserting (17)-(18) gives Continuity of saltwater requires the balance equation Combining (23) and (24) gives the evolution equation for the interface Equation (25) is solved subject to an initial condition where ζ 0 :R → [0, H ] satisfies certain smoothness conditions. In the absence of the integral term in (17), i.e., the Dupuit case, the term B(x; ζ ) drops from Eq. (25) and one is left with the nonlinear, doubly degenerate parabolic equation.
This equation received quite some attention in the mathematical literature. Its solutions behave as indicated in Fig. 1. Corresponding to an initial distribution as sketched in Fig. 1a, two free boundaries emerge in the (x, t)-plane, see Fig. 1b. One on the left, x = S 1 (t), separating the regions with only freshwater and fresh/saltwater. And one on the right, x = S 2 (t), separating the regions with fresh/saltwater and saltwater only. Existence, uniqueness, asymptotic behavior for large time and regularity and behavior of the free boundaries have been studied, see, for instance, Duijn and Hilhorst (1987), Van Duijn and Zhang (1988) and Bertsch et al. (1992).
A second main result of this paper concerns the validity of the Dupuit approximation, and consequently of Eq. (28). For this, we estimate the integral term in (17) (or in (42)) for q x and in (43) for q z . The results are given in Sect. 5.

Recovery of the Dupuit Discharge
The aim of this section is to derive (17) for q x and a similar expression for q z , both from (16). We start the derivation by assuming that for each t > 0, the interface ζ(x, t) is continuous in x with a piecewise continuous derivative, i.e., ζ(x, t) may have kinks.
Let (x p , z p ) be a fixed point in Ω which is either in freshwater (z p > ζ(x p , t)) or in saltwater (z p < ζ(x p , t)). At this point (x p , z p ), we want to construct a new expression for the specific discharge which clarifies the role and structure of the interface more explicitly.
For computational reasons, we introduce a small positive constant β, which we later send toward zero, such that ζ(x, t) = z p for x p − β < x < x p + β, see Fig. 2. By continuity  (16) the integration stretches from −∞ to +∞, we now split the integral into two parts: one where x s runs from −∞ to x p − β and one where x s runs from x p + β to +∞. Thus instead of (16) we consider the sum iΓ 2H where I 1 and I 2 denote the first and second part of the term between brackets in the integrand in (16). Clearly, the sum (29) tends to (16) as β → 0. We treat in detail the integration by parts of the term I 1 . The other term I 2 can be dealt with in a similar way. First we write the integrals containing I 1 in the form iΓ 2H where we use the notation ∂ζ . Because in the first term x s − x p < −β and in the second term x p − x s < −β, both the integrands can be considered as the limit of a uniformly convergent series. In fact, we have iΓ 2H Both integrals in (31) can be integrated by parts. The first integral in (31) yields Using expression (32) becomes To apply the same technique to the second integral in (31), we integrate the first term of the series separately, i.e., the term containing n = 0. Then we integrate by parts and substitute (33). The result is Next, we let β → 0 in (34) and (35). Using the identity ln(1 − e ±iα ) = ln 2 sin 1 2 α ± ig(α), we find, when adding (34) and (35), where 123 Here ζ x | x ± p = lim β→0 ζ x (x p ± β) and PV ∞ −∞ . . . dx s denotes the principal value of the integral at x p . In (37), and in other expressions to come, we use the plus-sign (+H ) when (x p , z p ) is chosen in the freshwater (z p > ζ(x p , t)) and we use the minus-sign (−H ) when (x p , z p ) is chosen in the saltwater (z p < ζ(x p , t)).
In a similar fashion, we rewrite the integrals containing the term I 2 in (29). The result is (skipping all details) The desired expression for the discharge is obtained after subtracting (37) and (39), i.e., Equation (40) is a general expression for the components q x and q z of the specific discharge at an arbitrary point (x p , z p ) in the aquifer. From this equation, a number of important conclusions can be drawn. First of all, it gives almost directly the shear flow at the interface. This follows from the choice of the sign in (37), as we shall see below. Secondly, it gives explicitly the influence of a jump in the derivative ζ x at x = x p . Having such a discontinuity, expression (37) shows that at the interface (z p → ζ(x p , t)), the components of the specific discharge become infinite, as if the fluids want to smooth out this discontinuity. Thirdly, it shows that the influence of the expressions which are related to the curvature of the interface, is accounted for by multiplying these terms by the weight function G and by integrating the result along the bottom of the aquifer. At a point ( which implies for the discharges (dropping the subscript p) and where Int(x, z; ζ ) is given by Expression (42) is the desired expression that we discussed in Sect. 2, see (17) and (18). The pointwise or local terms are the horizontal Dupuit discharges. Likewise expression (36) is written as where the vertical Dupuit discharge is given by While div q = 0 in Ω, this does not hold for q D . A straightforward differentiation gives Note that the Dupuit discharges do give the correct shear flow at the interface [see also (21)]. Letting z p → ζ(x p ) from above and from below yields the same integral contribution in (37). Hence, at z = ζ(x) we find and At points where ζ x = 0, for instance when only freshwater is present (ζ = 0) or when only saltwater is present (ζ = H ), the local Dupuit terms vanish from (42) and (43). At such points, the discharges are given by the integral terms only, which is in agreement with the findings of De Josselin de Jong (1981).

Validity of the Dupuit Discharge
Yortsos (1995) wrote a fundamental paper on the validity of the Dupuit approximation by considering the scaling: where L is a characteristic horizontal length scale. Introducing as a small parameter, he was able to recover the Dupuit approximation by applying an asymptotic expansion in terms of ε 2 . In his paper, Yortsos considered more general cases as we do here.
In expressions (42) and (43) the factor 1+ζ 2 x appears in the denominator. Applying scaling rules (50) and (51), with gives 1 + ε 2 u 2 x as denominator. Hence, in the expansion in terms of ε (or ε 2 ), this factor is removed from the expressions. We apply scaling (50)-(52) in Sect. 6 when considering the case of flat interfaces (ε << 1). In this section, we apply a more subtle approach in which we keep the influence of the denominator. We do this by analyzing (44). We show below in detail how this works for its real part, i.e., for q x . To obtain the real parts of (44), we first note that which we write as where Hence We estimate the kernels Im{ω − } and Re{ω + } in several steps.
Step 1 First consider ω − . Using (38) in (55) gives Note that, except for the sign of (x s −x), ω − for x s < x and for x s > x are each other complex conjugate. The structure of ω − is and hence Similarly Re{ω + } = ln(|T | · |B|) = ln |T | + ln |B| In (58) and (59), we have and Step 2 Let x s < x. Recall thatz = π z H , etc... Let a = e (x s −x) . Clearly 0 < a < 1. Then the expressions for T and B can be written as T = 1 + a e i(π −z+ζ ) and B = 1 + a e i(π −z−ζ ) .
In the complex plane, T and B are points on the circle C, see Fig. 3, with radius a and center (1,0). We make use of this figure to estimate Im{ω − } and Re{ω + }. - Step 3 Bound for Im{ω − }. With reference to Fig. 3, we observe that the maximum of |arg T | (or |arg B|) is attained when T (or B) is at the points K or M. The same holds for x s > x when a = e −(x s −x) . Hence for all x, x s ∈ R we have where a = e −|x s −x| .
Step 4 Bound for Re{ω + }. Again we use Fig. 3. Clearly Using (59) we obtain for all x, x s ∈ R where again a = e −|x s −x| . Note that estimates (63) and (64) do not depend on the position of the points T and B on the circle C. Therefore, they are uniform in z ∈ [0, H ] and also uniform with respect to the structure of ζ . With a = e −|x s −x| , let and Then we have obtained the estimate, recall (56) and dropped from now on the PV in the notation, for all x ∈ R and 0 ≤ z ≤ H . The kernels f and g are integrable in R since for each where c is a positive constant (0 < c < 2 ln 2 + 1). Furthermore, both kernels have the same asymptotic behavior for large |x s − x|. Let |x s − x| ≥ H . Then a ≤ e −π ≈ 1 20 . For such small a, we have and we estimate where we used (ii) and the monotonicity of f (x s − x 0 ). Since f (−1) = f (1) we further have Using +∞ −∞ f (x s − x) dx s ≤ 2 (apply(68) after scaling), (70) we find Similarly we find for I g Hence provided K is not too large. If K is large we need to extend the interval in (ii) of the Claim: say. This would give the factor e −2π in estimate (79). The right-hand side of (79) must be compared to the Dupuit discharge q D x from (18) after scaling (72) which is assumed O(1). This concludes the demonstration. For the vertical discharge component q z , one can proceed as above. Similar estimates are obtained and the Claim holds for q z as well. Details are omitted.
Example To clarify the conditions of the claim, we construct the following example. Let 0 < μ < 1 2 , α > 1 and drop the time dependence from the interface height : u = u(x). Then consider and see Fig. 4.  α, α). In this example, the conditions of the Claim, and hence the validity of Dupuit, hold for each point in the interval (α − 1, −α + 1). In the larger interval (−α, α) we have ∂ 2 u ∂ x 2 = 0 and This example also shows that for the Dupuit approximation to hold, it is not necessary to have very flat interfaces. If a = 3 2 , Dupuit holds in the interval (− 1 2 , 1 2 ) where du dx = 2 3 ( 1 2 − μ). It is therefore important to keep the factor 1 + ζ 2 x in the denominator of the Dupuit discharges (18), as suggested by de Josselin de Jong. In the next section, we consider the case of flat interfaces.

Flat Interfaces
We return to the scaling (50)-(52) introduced by Yortsos (1995). The assumption ε << 1 reflects indeed flat interfaces, since In the scaled variables, the Dupuit discharge (2.15) becomes . We show below that in the right-hand side of inequality (67), the f -term is O(ε 2 ) and the g-term O(ε 3 ). Applying (50)-(52) to (67) gives for x ∈ R and 0 ≤ z ≤ 1 123 where f ε is given by (65) and g ε by (66), but where now Since we find from (85) 2π In the integral, we set This gives and Heref is again given by (65) andg by (66), but now with Assuming that the integrals in (89) and (90) are bounded, uniformly in ε > 0, for (x, t) in a subdomain E of the (x, t)-plane, we have and consequently as ε ↓ 0. Consistent with the order of approximation, we may simplify q D x in (84) and in (93) by disregarding the term ε 2 u 2 x in the denominator. This gives The boundedness of the integrals in (89) and (90) where . Solutions of this equation should approximately describe the motion of the true interface, in particular for flat interfaces (ε << 1). Equation (95) has a particular solution, see Van Duijn and Zhang (1988) and De Josselin de Jong (1981), describing the redistribution of fresh and salt groundwater as a rotating line in space. It is given by the similarity solution  Using the definition (55) of ω + , we find Using (106) in expressions (42) and (43) gives the exact discharges in the different regions. Fig. 6, which is taken from De Josselin de Jong (1981), shows the discharges as well as the corresponding Dupuit discharges. Note that for x < S 1 (where ζ(x) = 0 and dζ dx = 0) and for x < S 2 (where ζ(x) = H and dζ dx = 0), the Dupuit discharges vanish. There the only flow contribution comes from (106a) and (106c).

Discussion and Conclusions
In this paper, we revisit the work of de Josselin de Jong on fresh-salt groundwater. In particular, we take his paper De Josselin de Jong (1981) as starting point for our analysis. In that paper, he considered fresh and salt groundwater in a horizontally extended aquifer. The scale of the problem is such that a sharp interface exists between the fluids. Based on earlier work De Josselin de Jong (1977), he constructs an expression for the horizontal and ver-tical discharge in terms of an integral along the fresh-salt interface. In doing so, he used the stream function formulation form De Josselin de Jong (1977). In this formulation, the stream function (and hence the discharge field) can be obtained for a given interface shape. In De Josselin de Jong (1981), he also considers a Dupuit-Forchheimer type approximation in which it is assumed that the horizontal discharge q x = q x (x, z) does not depend on z in each fluid. He verified this assumption using the stream function for a particular (piecewise linear) interface.
In this paper, we start from (16) as starting point. Applying integration by parts to (15), we recover the Dupuit-Forchheimer discharges as found by De Josselin de Jong (1981), see expressions (42)-(44). Moreover, we are able to give precise mathematical conditions for the validity of approximations (18) and (46). These conditions read as follows: If for some t > 0 and x 0 ∈ R we have i) ii) then the Dupuit approximation holds at x = x 0 , with an error that is small and that can be estimated. In other words, if (107) and (108) hold, the integral terms in (42) and (43) are small compared to the local Dupuit discharges (18) and (46). Note that (107) and (108) do not necessarily require flat interfaces. Therefore, the term (∂ζ /∂ x) 2 in the denominator in (18) and (46) cannot always be disregarded. In Sect. 5, we use the formulation from Sect. 3 and 4 to consider the case of flat interfaces. By flat we mean that ∂ζ /∂x = O(ε) and ε << 1. We give conditions for which we can estimate the order of convergence in terms of ε.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Derivation of Expression (15)
We obtain expression (15) as the solution of boundary value problem (12)-(14) in several steps. These steps are outlined here for completeness.
1. The stream function in the whole space R 2 due to a vortex of strength m(x s ): = Γ ∂ζ ∂ x (x s ) at the points (x s , ζ(x s )) is given by where r = (x − x s ) 2 + (z − ζ(x s )) 2 .