1 Introduction

In Sect. 2 of this paper, we solve the following harmonic problem for the complex potential \(w_1 ^{*}(z)=\phi _1 ^{*}+\hbox {i}\psi _1 ^{*}\) (holomorphic function in a complex plane \(z=x+\hbox {i}y\) with a cut):

$$\begin{aligned} \Delta \phi _1 ^{*} = 0,\quad \hbox { }\Delta \psi _1^*= 0,\quad - \infty<x<\infty ,\quad - \infty<y<\infty , \end{aligned}$$
$$\begin{aligned} \frac{\hbox {d} \phi _1 }{\hbox {d}x^*}=E^{*}\cdot \psi _1 \quad \hbox {for}\, y = 0, \quad - l<x<l, \end{aligned}$$
$$\begin{aligned} \psi _1 = 0 \quad \hbox {for}\,y = 0, \quad l<|x|, \end{aligned}$$
$$\begin{aligned} \frac{\hbox {d}w_1 }{\hbox {d}z^*}\rightarrow v_{\infty x} \quad \hbox { at }z\rightarrow \infty , \end{aligned}$$

where \(E^{*},l,\hbox { and }v_{\infty x}\) are given constants, and the linear relation between the stream function and the derivative of the velocity potential is called the Robin boundary condition. Similar boundary value problems have been solved in the theory of elasticity and aerodynamics (see, e.g. Golubev [1]). In Sect. 3, we solve an inverse problem for which the “cut” \({y=0, {\vert }x{\vert }<l}\), albeit remaining “thin”, has a variable aperture.

Physically, we consider a steady, saturated, one-phase, Darcian flow in a porous two-component composite made of a homogeneous isotropic rock (formation, aquifer) and a fracture perfectly spanning from the bedrock to caprock (Fig. 1) such that in any horizontal cross-section, the fracture contour \(\hbox {ABCB}_{\mathrm{b}}\hbox {A}\) is the same closed curve. The straight segment of the fracture axis AOC coincides with the abscissa axis \(\hbox {O}x\) of a Cartesian coordinate system. The fracture length is 2l. The aperture, d(s), where \(s (0<s<2l)\) is the arc length, satisfies the inequality \(d<<l\), i.e. the fracture is thin. Flow in each horizontal cross-section \(x^{*}y^{*}\hbox {O}\) of Fig. 1 is 2-D that is a limitation on the practical applicability of the obtained results to flows near natural and generated fractures, which are generally 3-D. However, the 2-D potential models are also widely used (see, e.g. van Harmelen and Weijermars [2]).

Applications of the flows considered below are as following:

  • Reservoir engineering (see, e.g. Bažant et al. [3]). The fracture is filled with a proppant, i.e. an artificial porous material injected into the fracture for keeping it open against geostatic pressure of thousands psi (e.g. in the Gulf hydrofracking practice, a dune sand is mixed with manufactured ceramics to make suitable proppants such that mixtures in uncompressed conditions have permeabilities of hundreds of darcies)

  • Environmental engineering (Klammler et al. [4]). The contour \(\hbox {ABCB}_{\mathrm{b}}\hbox {A}\) in Fig. 1 models a reactive barrier placed in aquifers; iron pellets are used as filling of the lens-shaped barriers; these barriers are not necessarily thin.

  • Subsurface hydrology. Porous lenses in rocks and soils are formed when loose sediments are accumulated by colmation, lessivage, suffusion or other hydrogeomorphological processes; \(\hbox {ABCB}_{\mathrm{b}}\hbox {A}\) is then a contour in a vertical cross-section, and the incident velocity is commonly vertical infiltration from above.

Refraction of Darcian flows on lenticular bodies, like one in Fig. 1, approximates them as ellipses or contours with corners in 2-D and as ellipsoids in 3-D Darcian flows (Adachi et al. [5], Chen [6], Dagan et al. [7], Kacimov and Obnosov [8], Mityushev and Adler [9], Lehr [10], Strack [11], Zhao et al. [12]). The refracted flows have drastically different kinematic characteristics, (pathlines-streaklines of marked particles, residence time in particular compartments of the subsurface porous composite, dispersion characteristics in associated transport of solutes, etc., see, e.g. Kacimov and Obnosov [13, 14]) compared with “effectively” homogenized media. This has dramatic consequences in applications: for example, in natural or managed aquifer recharge of unconfined aquifers in desert zones, much less water arrives at the water table from the infiltrated ground surface than is expected according to hypothesized homogeneous or simply stratified vadose zone models. The imbedded lenses increase the length of the pore water path, which occurs to be much more “winding” compared with “standard” stratification patterns conceptualized by hydrogeologists and vadose zone hydrologists. Consequently, the travel time of soil water increases such that plant roots manage to intercept it and transpire back to the atmosphere.

Fig. 1
figure 1

Formation with a perfect vertical fracture

Full size image

The hydraulic conductivities of the ambient rock and of fracture’s filling in Fig. 1 are \(k_{b}\) and \(k_{f}\), correspondingly. We focus on the case \(k_{b}\ll k_{f}\). If the fracture is empty, then \(k_{f}\) is proportional to \(d^{3}\) (Pilatovskii [15], Zhang et al. [16]). We study planar fractures having two axes of symmetry and select the origin of a Cartesian coordinate system \(x^{*}Oy^{*}\) coinciding with the centre O of the fracture (Fig. 1).

Far away from the fracture, a unidirectional ambient flow is arbitrarily oriented with respect to the fracture. The Darcian velocity vector at infinity (point D in Fig. 1), \({\vec {V}}^{*}_\infty \), has two components \((v^{*}_{\infty x} , v^{*}_{\infty y})\) along the \(x^{*}\)- and \(y^{*}\)-axis, respectively. This “far-field” condition corresponds to a mathematical dipole at infinity (Strack [11]). In the vicinity of the fracture (“near-field”), flow is essentially 2-D. A part of the incident flow is intercepted by the fracture, fluid moves along its axis and seeps out, back to the ambient rock. Two separatrices, each consisting of two branches (\(\hbox {D}_{1}\hbox {S}_{\mathrm{u}}-\hbox {S}_{\mathrm{u}}\hbox {D}_{2}\) and \(\hbox {D}_{3}\hbox {S}_{\mathrm{d}}-\hbox {S}_{\mathrm{d}}\hbox {D}_{4}\) in Fig. 1) are formed, as investigated by Kacimov et al. [17], Kacimov and Obnosov [8] for applications to problems of full refraction on elliptical lenses. The width of the capture zone is 2*in. In this paper, we study two fracture geometries: \(d(s)=\delta =const\) and fractures of variable apertures, for which d(s) is obtained from solution of an inverse problem (with a given flow rate along fracture’s axis). The main accomplishment is in an analytical treatment of the Robin boundary condition along the contours of these fractures and in illustration of applicability-extension of the methods developed in mechanics of ideal fluids (thin airfoils) to mathematically equivalent problems in subsurface mechanics.

2 Fractures of constant aperture

In this section, we consider a fracture of a small constant aperture \(\delta \). We assume that the fracture is a mathematical slot (cut). The value of \(k_{f}\) is so high and the rock is so tight that the dimensional quantity \(E^*=k_b /(k_f \cdot \delta )=\hbox {const}\) is relatively small even for thin fractures.

2.1 Formulation of a boundary-value problem

Complex variables are used in modelling of steady Darcian potential flows near fractures (see, e.g. Liolios and Exadaktylos [18], Strack [11], Thiele et al. [19], van Harmelen and Weijermars [2], Weijermars and Alves [20]).

We introduce a complex physical coordinate \(z^*=x^*+\hbox {i}\cdot y^*\), the longitudinal axis of the fracture coincides with the real axis, \(y^{*}=0\). A complex potential is defined as \(w^{*}=\phi ^{*}+\hbox {i}\psi ^{*}\), where \(\varphi ^{*}\) is the velocity potential, i.e. the Darcian velocity is \({\vec {V^{*}}}(u^{*},v^{*})=\nabla \varphi ^{*}\) and \(\psi ^{*}\) is the stream function which is complex-conjugated with \(\varphi ^{*}\) according to the Cauchy–Riemann conditions (see, e.g. Strack [11]), \(u^{*}\hbox { and }v^{*}\) are the components of the velocity vector along \(x^{*}\) and \(y^{*}\) axes.

We use the principle of superposition (see, van Harmelen and Weijermars [2]), i.e. decouple \(w^*(z^*)\) as \(w^*(z^*)=w_1 (z^*)+w_2 (z^*)\), where \(w_2 (z^*)=-\hbox {i}\cdot v^{*}_{\infty y} \cdot z^*\) and therefore \(w_1 (z^*)=w^*(z^*)-w_2 (z^*)\). For \(w_1 (z^*)\), the velocity at infinity \({\vec {V}}^{*}_\infty \) has two components \((v^{*}_{\infty x} , 0)\), i.e. this vector is oriented along the real axis of Fig. 2a. For \(w_2 (z^*)\), the incident velocity far from the fracture is oriented along the imaginary axis. Obviously, both \(w_1 (z^*)\hbox { and }w_2 (z^*)\) are analytic functions.

Fig. 2
figure 2

A left upper quarter of the a flow domain, b complex potential domain and c a reference plain

Full size image

We have to find \(w_1 (z^*)\) in the plane \(z^*\) with a cut \(-\,l<x^{*}<l\), \(y^{*}=0\) such that \({\hbox {d}w_1}/{\hbox {d}z^{*}}\rightarrow v_{\infty x} \) at \(z^{*}\rightarrow \infty \). Analogously, for the case of an incident flow at infinity parallel to the ordinate axis, for \(w_2 (z^{*})\) the asymptotics is \({\hbox {d}w_2 }/{\hbox {d}z^{*}}=-\hbox {i}\cdot v_{\infty y}\) at \(z^{*}\rightarrow \infty \).

The complex potential is scaled up to an arbitrary constant and due to symmetry we select \(\upvarphi _1 =0\) along the imaginary axis \(x^{*}=0\) in Fig. 2a. Also, we set \(\uppsi _1 =0\) along \(y^{*}=0\), \(x^{*}<-l, x^{*}>l\). Along the fracture contour \({BCB}_{r}\), \(-l<x^{*}<l\), \(y^{*}=0\) the condition holds:

$$\begin{aligned} \frac{\hbox {d} \phi _1 }{\hbox {d}x^{*}}=E^{*}\cdot \psi _1 . \end{aligned}$$
(1)

The Robin boundary condition (1) was first derived for Darcian flows near thin fractures in [15] . Physically, Eq. (1) implies that within the fracture the “internal” flow is 1-D along the fracture axis, but the flow rate of this in-fracture Darcian flow increases from zero at the left tip A (Fig. 2a) to a maximal value at cross-section OB and then decreases again to zero at the left tip C (Fig. 1). Correspondingly, the complex potential of the “external” flow in the aquifer (formation), i.e. in the whole z-plane of Fig. 2a is conjugated with the “internal” flow. Although the aperture of the fracture in Fig. 2a is physically finite (and in Sect. 3 is variable with the x-coordinate), in the mathematical model a zero-aperture is assumed and the “external” flow problem is solved in the whole plane z with a “contact” boundary condition along the cut \({{\vert }x{\vert }<l}\) in the vernacular of the theory of elasticity.

This condition used for “thin objects”, either highly- or low-permeable (fractures and cutoff walls) makes the problem complicated for analytical treatment, with only few solutions known in the literature (Anderson [21], Kacimov and Al-Maktoumi [22], Kacimov and Obnosov [23, 24], Strack [11], Van Der Veer [25], Yakimov and Kacimov [26]).

We introduce dimensionless variables \((z, x, y)=(z^{*}, x^{*}, y^{*})/l\), \((w, \varphi , \psi )=(w^{*}_1 , \varphi ^{*}_1, \psi ^{*}_1 )/(v^{*}_{\infty x} \cdot l),(V,u,v)=(V^{*},u^{*},v^{*})/v^{*}_{\infty x}\), \(E=E^{*}\cdot l\), such that the parameter

$$\begin{aligned} E=k_b \cdot l/(k_f \cdot \updelta ) \end{aligned}$$
(2)

varies from 0 (the limit of an ideally conducting fracture of zero hydraulic resistance, the case studied by van Harmelen and Weijermars [2]) to infinity (an impermeable wall, e.g. a geological vein cemented by precipitation of carbonates from seeping groundwater).

Thus, we have the following boundary value problem (BVP). Determine the complex potential \(w(x)=\phi (x, y)+\hbox {i}\cdot \psi (x, y)\), which is analytic in the whole plane z and satisfies the conditions along the cut:

$$\begin{aligned} \frac{\partial \upvarphi }{\partial x}=E\cdot \uppsi \quad \hbox { at }-1<x<1,\quad y=0, \end{aligned}$$
(3a)

and a symmetry condition:

$$\begin{aligned} \uppsi =0\quad \hbox { at }y=0,\quad |x|>1. \end{aligned}$$
(3b)

The behavior of w(z) at infinity is \({\mathop {\lim }\nolimits _{z\rightarrow \infty }} ({\hbox {d}w}/{\hbox {d}z})\rightarrow 1\). Figure 2b shows the complex potential domain \(\hbox {G}_{w}\) which corresponds to the upper left quarter of the physical domain \(\hbox {G}_{\mathrm{z}}\) in Fig. 1a. In \(\hbox {G}_{\mathrm{w}}\), the image of the quarter of the contour AB of the fracture is an unknown curve. For an ambient flow oriented along the longitudinal axis of the fracture (Fig. 2a), the stagnation point \({S}_{\mathrm{u}}\) in Fig. 1 coincides with point B and the branch of the separatrix \({D}_{1}{S}_{\mathrm{u}}\) in Fig. 1 is imaged by a ray \({D}_{1}{B}\) shown by a dashed line in Fig. 2b.

We introduce an auxiliary complex variable \(\zeta =\xi +\hbox {i}\cdot \eta =r\cdot \hbox {e}^{\hbox {i}\cdot \theta }\) and map conformally the exterior of the cut \(\hbox {ABCB}_{\mathrm{b}}\hbox {A}\) (\(-1<x<1\), \(y=0\)) in Fig. 2a on the exterior of the circle \(|\upzeta |>1\) by the Zhukovsky function (Fig. 2c):

$$\begin{aligned} z=\frac{1}{2}\left( {\upzeta +\frac{1}{\upzeta }} \right) , \end{aligned}$$
(4)

which is inverted by

$$\begin{aligned} \upzeta =z+\sqrt{z^{2}-1} \end{aligned}$$
(5)

(see, e.g. [1]).

In the reference plane \(\upzeta \), the function \(w(\upzeta )\) is symmetric, such that for \(r\ge 1\), we have the following conditions:

$$\begin{aligned}&\hbox {DA and CD:}\quad \uppsi =0\hbox { along the real axis }\upeta =0,\nonumber \\&\hbox {BD:}\quad \upvarphi =0\hbox { along the imaginary axis }\upxi =0. \end{aligned}$$
(6)

From \({\hbox {d}w}/{\hbox {d}z}\rightarrow 1\) and Eq. (4), the far-field condition is

$$\begin{aligned} \frac{\hbox {d}w(\zeta )}{\hbox {d}\zeta }\rightarrow \frac{1}{2}\hbox { at }\upzeta \rightarrow \infty . \end{aligned}$$
(7)

The boundary \(\hbox {ABCB}_{\mathrm{b}}\hbox {A}\) of the z-plane is imaged by the circumference \(r=1\) in the reference \(\upzeta \)-plane, with the following correspondence of variables in the two planes:

$$\begin{aligned} x=\cos \uptheta . \end{aligned}$$
(8)

Therefore, using Eq. (8) the condition (3a) is reduced to

$$\begin{aligned} \frac{\hbox {d}\phi (\theta )}{\hbox {d}\theta }=-E\cdot \psi (\theta )\cdot \sin \theta . \end{aligned}$$
(9)

2.2 Robin’s boundary condition transformed to an infinite system of linear equations

We search for \(w(\upzeta )\) as a series, in which the symmetry conditions (6) make the expansion coefficients \(A_{i}\) real numbers. The expansion has only terms with odd exponents of \(\zeta \):

$$\begin{aligned} w(\upzeta )=A_0 \cdot \upzeta +\frac{A_1 }{\upzeta }+\frac{A_2 }{\upzeta ^{3}}+\cdots +\frac{A_n }{\upzeta ^{2n-1}}+\cdots \,\,. \end{aligned}$$
(10)

Obviously, \(A_0 =1/2\) owing to Eq. (7).

From Eq. (10), for \(\upvarphi (\uptheta )\) and \(\uppsi (\uptheta )\) along the circle ABC in the \(\zeta \)-plane, we have

$$\begin{aligned} \upvarphi (\uptheta )= & {} (A_0 +A_1 )\cdot \cos \uptheta +A_2 \cdot \cos (3 \uptheta )+\cdots +A_n \cdot \cos [(2n-1 )\cdot \uptheta ]+\cdots \,, \end{aligned}$$
(11)
$$\begin{aligned} \uppsi (\uptheta )= & {} (A_0 -A_1 )\cdot \sin \uptheta -A_2 \cdot \sin (3 \uptheta )-\ldots -A_n \cdot \sin [(2n-1 )\cdot \uptheta ]-\ldots \,. \end{aligned}$$
(12)

Thus, we have to find \(A_{n}\) in Eqs. (11)–(12), which satisfy Eq. (9).

To the right-hand side of Eq. (9), i.e. to \(\uppsi (\uptheta )\cdot \sin \uptheta \), we apply the trigonometrical formula:

$$\begin{aligned} \sin (m\uptheta )\cdot \sin \uptheta =\frac{1}{2}\Big \{\cos [(m-1)\uptheta ]-\cos [(m+1)\uptheta ]\Big \},\hbox { or }\nonumber \\ \sin [(2n-1)\uptheta ]\cdot \sin \uptheta =\frac{1}{2}\Big \{\cos [(2n-2)\uptheta ]-\cos (2n\uptheta )\Big \}, \end{aligned}$$

such that

$$\begin{aligned} \uppsi (\uptheta )\cdot \sin \uptheta= & {} \frac{1}{2}\Big \{-(A_1 -A_0 )+(A_1 -A_0 -A_2 )\cos (2\uptheta )+(A_2 -A_3 )\cos (4\uptheta )+\cdots \nonumber \\&\cdots +(A_n -A_{n+1} )\cos (2n\uptheta )+\cdots \Big \}. \end{aligned}$$
(13)

The left-hand side of Eq. (9) is represented as

$$\begin{aligned} \frac{\hbox {d}\phi (\theta )}{\hbox {d}\theta }=-\Big \{(A_1 +A_0 )\sin \theta +3A_2 \sin (3 \theta )+\ldots +(2n-1)A_n \sin [(2n-1 )\theta ]+\ldots \Big \}. \end{aligned}$$
(14)

Next, we put the expansions (13) and (14) into Eq. (9), multiply both sides by \(\sin (m\uptheta )\) and integrate from 0 to \(\uppi \). We use the following orthogonality conditions:

$$\begin{aligned}&\int _0^\pi {\sin (k\theta )\cdot \sin (m\theta ) \hbox {d}\theta } = \frac{1}{2}\int _0^\pi {\Big \{\cos [(k-m)\theta ]-\cos [(k+m)\theta ]\Big \} \hbox {d}\theta } =\left\{ {{\begin{array}{ll} \pi /2 &{}\hbox { for } k=m \\ 0 &{}\hbox { for } k\ne m \\ \end{array} }} \right. \end{aligned}$$
(15)
$$\begin{aligned}&\int _0^\pi {\cos (2n\theta )\cdot \sin (m\theta ) \hbox {d}\theta } = \frac{1}{2}\int _0^\pi {\Big \{-\sin [(2n-m)\theta ]+\sin [(2n+m)\theta ]\Big \}\hbox {d}\theta } =\frac{2m}{m^{2}-(2n)^{2}} \,\hbox { for odd }\, m. \end{aligned}$$
(16)

Consequently, for odd m, i.e. for \(m=2n-1, n=1, 2, \ldots \), Eq. (9) gives

$$\begin{aligned} \int _0^\pi {\frac{\hbox {d}\phi (\theta )}{\hbox {d}\theta }\cdot \sin (m\theta ) \hbox {d}\theta } =-\int _0^\pi {(E\cdot \psi (\theta )\cdot \sin \theta )\cdot \sin (m\theta ) \hbox {d}\theta } , \end{aligned}$$

or upon putting into this relation Eqs. (13) and (14) along with Eqs. (15) and (16), we get an infinite system of equations

$$\begin{aligned} -\tilde{B}_n \cdot m\cdot \frac{\pi }{2}=-E\cdot \sum _{j=0}^\infty {(B_j -B_{j+1} )\frac{m}{m^{2}-(2j)^{2}} },\quad m=2n-1,\quad n=1, 2, 3, \ldots , \end{aligned}$$
(17)

where we use the notations

(18)

It is noteworthy that series expansions like ones above have been widely used in solving BVPs to harmonic functions in canonic domains (see, e.g. Tikhonov and Samarskii [27]) but mainly for Dirichlet and/or Newman, rather than Robin’s boundary condition.

2.3 Proof of regularity of infinite system

Systems of an infinite number of linear equations, like Eq. (17), with an infinite number of unknowns emerged in several applications in mechanics of fluids and elastic bodies, with a theoretical compendium in Kantorovich and Krylov [28] . Aravin and Numerov [29], Kasimov and Tartakovsky [30] applied this theory to Darcian flows.

System (17) can be re-written in the form [28] (see Chapter I, §2, p. 30 in their book)

$$\begin{aligned} B_n =\frac{2E}{\uppi }\sum _{k=1}^\infty {c_{nk}^*\cdot B_k } +b_n , \quad n=1, 2, 3, \ldots , \end{aligned}$$
(19)

\(b_1 =-1, \quad b_n =0\) at \(n>1\), or

$$\begin{aligned} A_n =\frac{2E}{\uppi }\cdot \left( {\sum _{k=1}^\infty {c_{nk}^*\cdot A_k } } \right) +b_n , \quad n=1, 2, 3, \ldots , \end{aligned}$$
(20)

where

$$\begin{aligned} b_1 =-\frac{1}{2}+\frac{4E}{3\uppi },\quad b_n =0\quad \mathrm{at}\; n>1, \end{aligned}$$
$$\begin{aligned} c_{nk}^*=\frac{1}{(2n-1)^{2}-4k^{2}}-\frac{1}{(2n-1)^{2}-4(k-1)^{2}}=\frac{4(2k-1)}{[(2n-1)^{2}-4(k-1)^{2}]\cdot [(2n-1)^{2}-4k^{2}]}. \end{aligned}$$

Using Wolfram’s [31] Mathematica routine Sum, we get \(\sum \nolimits _{k=1}^\infty {c_{nk}^*} =-{{1}/{(2n-1)^{2}}}\). It is noteworthy that \(c_{nn}^*<0\), and \(c_{nk}^*>0\) at \(n\ne k\). Therefore,

$$\begin{aligned} \sum \limits _{k=1}^\infty {\left| {c_{nk}^*} \right| } =\left( {\sum \limits _{k=1}^\infty {c_{nk}^*} } \right) -2c_{nn}^*=\frac{64n^{3}-112n^{2}+64n-11}{(2n-1)^{2}(4n-3)(4n-1)}. \end{aligned}$$

Obviously, the maximum

$$\begin{aligned} \underbrace{\max }_n\left( {\sum \limits _{k=1}^\infty {\left| {c_{nk}^*} \right| } } \right) =\sum \limits _{k=1}^\infty {\left| {c_{1k}^*} \right| } =\frac{5}{3}. \end{aligned}$$

In [28] (see their Eq. 25, Chapter I, §2, p. 37) entirely regular infinite systems of linear equations defined ones for which the sum of absolute values of the coefficients of the system is less than \(1-\upvarepsilon \), \(\upvarepsilon \) is a small positive constant.

At

$$\begin{aligned} E<\frac{\uppi }{2\cdot \max \left( {\sum \limits _{k=1}^\infty {\left| {c_{nk}^*} \right| } } \right) }=\frac{3\uppi }{10}\approx 0.942 \end{aligned}$$
(21)

our system (19) or (20) is entirely regular.

In [28] (see their Theorems Ia, IIa, and Comment on Theorem IIa, Chapter I, §2, pp. 38–39) proved that entirely regular systems have a unique bounded solution, which can be found by the method of iterations. Moreover, Theorem IVa of [28] guarantees that such solution can be found by the method of reduction, viz. the infinite system is truncated to N equations and N unknowns and at a sufficiently large N the solution of a reduced (finite) system practically coincides with the solution of the original infinite system, i.e. at \(N\rightarrow \infty \) the solution of a finite system converges.

We use this reduction approach, i.e. truncate the system and solve a finite system by iterations with a subsequent testing that at a sufficiently large N the solution is not sensitive to N. In other words, we chop Eq. (17) to:

$$\begin{aligned} \tilde{B}_n^{(k)} =\frac{2\cdot E}{\pi \cdot m}\cdot \sum _{j=0}^N {\left( B_j^{(k)} -B_{j+1}^{(k)}\right) \frac{m}{m^{2}-(2j)^{2}}},\quad m=2n-1,\quad n=1, 2, 3, \ldots , N. \end{aligned}$$
(22)

According to Eq. (18) the iterations are generated in the following manner:

$$\begin{aligned} B_1^{(k+1)} =\tilde{B}_1^{(k)} -2A_0 , \quad B_n^{(k+1)} =\tilde{B}_n^{(k)}\quad \hbox { for }\quad n= 2, 3,\ldots , N. \end{aligned}$$
(23)

where the superscript \(^{(k)}\) designates the k-th iteration. As the initial iteration, we select \(\tilde{B}_1^{(0)} =A_0 \), \(\tilde{B}_n^{(0)} =0\) at \( n= 2, 3,\ldots N\). Correspondingly, after solving the system, in Eqs. (10) or (11), (12) (remembering that \(w(z)=w(\upzeta (z))\), together with Eq. (5), we have

$$\begin{aligned} A_1 =B_1 +A_0 , \quad A_n =B_n\quad \hbox { for } n= 2, 3, \ldots , N. \end{aligned}$$
(24)

We recall that the whole complex potential at an arbitrary orientation of the incident flow, involves \(w_2\). For the latter, in the introduced dimensionless variables the “far field” condition is \(\phi _2 =v^{*}_{\infty y} \hbox { }y/(v^{*}_{\infty x} \cdot l), \psi _2 =-v^{*}_{\infty y} \hbox { }x/(v^{*}_{\infty x} \cdot l).\)

In order to speed up the iterations, we use the following algorithm (similar to one implemented in [17, 26]), which involves three consecutive iterations for each coefficient. Specifically, for \(B_n^{(k+1)}\), \(B_n^{(k-1)}\), \(B_n^{(k)}\), we use the formula:

$$\begin{aligned} B_n^*=\frac{\left( B_n^{(k)}\right) ^{2}-B_n^{(k+1)} \cdot B_n^{(k-1)} }{2B_n^{(k)} -B_n^{(k+1)} -B_n^{(k-1)} }, \quad n= 1, 2, 3, \ldots , N, \end{aligned}$$
(25)

Below we describe the results of numerical experiments with variation of a physical parameter E, and two mathematical parameters, viz. N, relative error of computed coefficients, \(e_{A}\), and accuracy of meeting the boundary condition (9), \(e_{9}.\) Numerical experiments showed that primitive iterations (23) converge at \(E\le 1.1\). We recall that according to the theory in [28] the iterations converge at \(E<0.942\) according to Eq. (21).

For instance, for \(E=1\) and selected \(e_{A} =10^{-10}\) evaluation of coefficients \(A_i \) required 140 iterations. In the boundary condition (9), the accuracy of computations of the derivatives depends on the number of truncated terms: for instance, for \(N=500\), the accuracy \(e_{9} = 10^{-8}-10^{-9}\). The algorithm of an accelerated convergence (22), gives the same accuracy quantified by a dyad (\(e_{A, }e_{9}\)) in 9–11 iterations for \(100<N<1000\). Moreover, the “three-point” iterations (22) converge at E values as high as 52. For example, at \(E=52\), the accuracy \(e_{A} =10^{-8}\) for \(A_i\) and \(e_{9} =10^{-6}-10^{-7}\) in Eq. (9) (retaining \(N=500\) in the truncated series) is attained with 34 iterations in Eq. (22).

Fig. 3
figure 3

a Flow net near rectangular fracture for \(v_\infty =0.0,E=0.6,N=500\). b Flow net for \(v_\infty =0.5,E=0.6,N=500\)

Full size image

The model in [15] assumes 1-D flow inside the fracture (along AOC in Fig. 2a) with accretion-losses from the ambient rock. This 1-D flow in the fracture is—in a sense—similar to the Dupuit-Forchheimer 1-D phreatic flow in groundwater hydrology (Strack [11]). For an ideally conducting (no resistance) fracture \(E=0\), AB in Fig. 2b is an equipotential contour \(\upvarphi =0\) (\(\uppsi \) varies from 0 to 1 along AB). For an “impermeable” fracture (barrier), \(E=\infty \), AB in Fig. 2b is a streamline \(\uppsi =0\) and \(\upvarphi \) varies from − 1 to 0 along AB. The limit of large E is tackled by a model of a low-permeable barrier (cutoff wall), the cross-flow through which is proportional to the difference of the piezometric head at opposite sides of the cut along which the Robin boundary condition is mathematically similar to Eq. (3a) ([21, 26]).

As an example, we select physical parameters corresponding to \(E=1\). Then at the two tips of the fracture, A and C, \(\varphi =\pm 0.42\) (correspondingly) and \(\psi _B =Q_{f} =0.52\) (see, Fig. 2b) at the middle point B. Therefore, the perturbations of the incident 1-D flow caused by this fracture are quite significant and the quantity \(Q_{f}\) of the fluid intercepted by the fracture is high. Correspondingly high is the width \({2*in}=2{Q}_{f}/k_{b}\) of the “capture zone” of the fracture. This width equals the distance between the branches of separatrices at infinity (Figs. 1 and 2a). In the case of a fracture having \(E=52\), the perturbation of the longitudinal ambient flow is much less: \(|\varphi _{A,C}| =0.962\) at the tips and \(\psi _B =0.019\). From this comparison of two values of E, we can infer the following assessment for hydrofracking. The 52-fold reduction of the fracture opening \(\updelta \), provided the other parameters determining E (i.e. \(l, k_{f}, k_{b})\) remain the same, results in a 27-fold decrease of the quantity of the fluid intercepted by the fracture in a unidirectional ambient flow.

Using the ContourPlot, Re and Im routines of Mathematica, in Fig. 3a, we plotted the flow net for the case of \(v_\infty =0.0,E=0.6,N=500\). The equipotential lines span from \(\varphi =-1.2\) (curve 1) to \(\varphi =1.2\) (curve 2) with an increment of 0.2. The streamlines are plotted from \(\psi =0\) (abscissa axis) to \(\psi =1.6\) (curve 3) with an increment of 0.2. The streamline \(\psi =Q_f \) delineates the part of the incident flow which is intercepted by the fracture. For this example, the separatrix is \(Q_f =0.65\) (curve 4). It is clear from Fig. 3a that the streamlines which enter and exit the fracture are not orthogonal to it, as would be in the case of a fracture of infinite conductivity; the slope of streamlines at the intersection with the fracture gradually decreases from the left tip of the fracture to its middle.

In Fig. 3b, we present the flow net for the case of \(v_\infty =0.5,E=0.6,N=500\). The equipotential lines span from \(\varphi =-1.1\) (curve 1) to \(\varphi =1.9\) (curve 2) with an increment of 0.1. The streamlines are plotted from \(\psi =-0.62\) (curve 3) to \(\psi =2.18\) (curve 4) with an increment of 0.2. It is noteworthy that, unlike the case of Fig. 3a, in b the two-branch streamlines indicating inflow-outflow into and from the fracture are asymmetric: for instance, the streamline \(\psi =0.78\) arrow-indicated as two branches of curve 5.

3 Pilatovskii’s “inverse” problem

Pilatovskii [15] proposed an inverse method for calculation of flow in porous formations with fractures having variable apertures and finite (non-zero) resistance to flow inside them. Later, the ideas of [15] have been extended in a series of papers by Abdurakhmanov [32], Abdurakhmanov and Alishaev [33]. Abdurakhmanov and Alibekov [34] applied the method of [15] to mathematically equivalent problems of heat conduction through solid bodies with fractures.

In the inverse approach applied to seepage problems, a hydraulic characteristic (distribution of a hydraulic gradient, velocity potential, stream function, etc.) of an object under design (earth dam, soil channel, drain, etc.) is specified, the geometry of the object is found as a part of solution of BVP (see, e.g. Ilyinsky and Kacimov [35], Ilyinsky et al. [36]).

Pilatovskii [15] specified the flow rate, \(\omega \,(x)\), inside the fracture along its axis and found \(\delta (x)\). Without any loss of generality, we assume that the external velocity in Fig. 2a is oriented along the abscissa axis and, as in Sect. 2, consider fractures symmetric with respect to the \(\hbox {O}x\) and \(\hbox {O}y\) axes. Then the flow rate function in [15] obeys the properties:

$$\begin{aligned} \omega ^{*}(-x^{*})=\omega ^{*}(x^{*}),\omega ^{*}(\pm l)=0. \end{aligned}$$
(26)

This arbitrary function belongs to the class \(\hbox {C}^{1}\), i.e. its derivative obeys the Holder condition that is sufficient for the existence of the Cauchy integrals below.

Obviously, for a practically occurring fractures with a single maximum, \(2\delta { }_0=2\delta (0)\), of the fracture opening at point O (Fig. 2a) the inequality \(\omega (x{ }_1)<\omega (x{ }_2)\) for \(|x{ }_1|>|x{ }_2|\) holds and the maximum, \(2Q_{f}\), of \(\omega (x)\) is attained in the middle of the fracture at \(x=0\).

In this section, dimensionless variables are introduced as: \((z, x, y,\delta )=(z^*, x^*, y^*,\delta ^*)/l\), \(V_0 =V_0^{*} /k_b ,(w, \varphi , \psi ,\omega )=(w^*_1 , \varphi ^*_1 , \psi ^*_1 ,\omega ^*)/(k_b \cdot l)\).

For a planar fracture obeying Eq. (1), Pilatovskii [15] derived an integro-differential equation:

$$\begin{aligned} \frac{\omega (x)}{\delta (x)}=V_0 +\frac{1}{2\pi }\frac{\hbox {d}}{\hbox {d}x}\int \limits _{-1}^1 {\frac{\omega (u)}{u-x}} \hbox {d}u,\hbox { }\quad -1\le x\le 1, \end{aligned}$$
(27)

which interrelates the dyad [\(\delta (x)\), \(\omega (x)\)], i.e. the geometric and hydraulic characteristics. The Cauchy integral in Eq. (27) is singular, and we evaluated it in the sense of V.P., either analytically or using the routine NIntegrate,

Method->{PrincipalValue} from [31].

It is noteworthy that Eq. (27) is modulo notations, a special case of the Prandtl equation in aerodynamics of thin airfoils. In the Prandtl flow of an ideal fluid past a wing of an aircraft, \(\omega (x)\) is a distributed singularity (circulation) along an airfoil chord. A fictitious fluid circulates inside the Prandtl airfoil. Physically, flow inside the airfoil (which is a mathematical separatrix or a physically impermeable flow boundary) is discarded, while in applications to Darcian flows, both the interior and exterior of separatrix-confined subdomains are physically meaningful analytic elements (see, e.g. Strack [11]).

Obviously, Eq. (27) can be re-written as

$$\begin{aligned} \delta (x)=\frac{\omega (x)}{V_0 -\frac{1}{2\pi }\int \nolimits _{-1}^1 {\frac{\omega ^{\prime }(u)}{u-x}} \hbox {d}u},\quad -1\le x\le 1. \end{aligned}$$
(28)

Owing to fractures’ symmetry, we consider the upper half-plane in Fig. 2a (\(y>0\)). Pilatovskii [15] investigated a special case of \(\omega _P (x)=Q_{fP}(1-x^2),{\omega _P}^{\prime }(x)=-2Q_{fP}x\) for which Eq. (28) for the contour of the fracture is reduced to:

$$\begin{aligned} \delta _P (x)=\frac{Q_{fP} (1-x^{2})}{V_0 +\frac{Q_{fP} }{\pi }\left[ {x\ln \frac{1+x}{1-x}-2} \right] },\quad \hbox { }\delta _{0P} =\frac{Q_{fP} }{V_0 -\frac{2Q_{fP} }{\pi }}. \end{aligned}$$
(29)

Below we consider an arbitrary \(\omega (x)\), which is written as the following series expansion of its derivative, i.e. the numerator of the integrand in Eq. (29):

$$\begin{aligned} \omega ^{\prime }(x)=\omega ^{\prime }_P (x)+Q{}_{fP}^{*}\sum _{n=1}^\infty {b{ }_{2n}} U_{2n} (x),U_n (x)=\sin [n\arccos (x)],\hbox { }-1\le x\le 1. \end{aligned}$$
(30)

Obviously, at \(b_{2n} \equiv 0\), we get Pilatovskii’s fracture.

The coefficients \(b_{2n}\) in the series involving Chebyshev’s functions are real constants and upon evaluation of the Cauchy integral in Eq. (28) via Chebyshev’s polynomials of the first kind, \(T_{n}(x)=\hbox {cos}[n\hbox { arccos}(x)]\) (see, e.g. Kacimov [37]), we obtain

$$\begin{aligned} \delta (x)=\frac{(1-x^{2})+\sum \nolimits _{n=1}^\infty {b_{2n} \frac{2n\sqrt{1-x^{2}}T_{2n} (x)-xU_{2n} (x)}{4n^{2}-1}} }{\frac{V_0 }{Q{ }_{fP}}+\frac{1}{\pi }\left[ {x\ln \frac{1+x}{1-x}-2} \right] +\frac{1}{2}\sum \nolimits _{n=1}^\infty {b{ }_{2n}} T_{2n} (x)},\quad \hbox { }-1\le x\le 1. \end{aligned}$$
(31)

If the shape of a fracture is fixed then \(b_{2n}\) can be found as for Prandtl’s thin wings . Namely, truncation of the series up to N, selection of N collocation points \(x_{2n}\) and solution of a system of N linear equations with respect to \(b_{2n}\) [obtained from Eq. (31) with the RHS equal to constants \(\delta (x_{2n})\)] can be carried out similarly to what we did in Sect. 2.

Following the inverse procedure of [15], we specify \(b_{2n}\) and use Eq. (31) for reconstruction of the geometry, \(\delta (x)\). For example, if we retain only one term \(b_{2}\) in the series then Eq. (31) yields

$$\begin{aligned} \delta _P (x)=\frac{(1-x^{2})-\frac{2}{3}b_2 (1-x^{2}){ }^{3/2}}{\frac{V_0 }{Q{ }_{fP}}+\frac{1}{\pi }\left[ {x\ln \frac{1+x}{1-x}-2} \right] +b_2 (x^{2}-0.5)}. \end{aligned}$$
(32)
Fig. 4
figure 4

Shapes of curvilinear half-fractures, Pilatovskii’s contour (curve 1, \(b_{2}=0\)), and its generalizations \(b_{2}=0.2\) and 0.8 (curves 2 and 3)

Full size image
Fig. 5
figure 5

Streamlines near a curvilinear fracture

Full size image

In Fig. 4, the curves 1–3 correspond to the fracture contours at \(b_{2}=0, 0.2, 0.8\) plotted by Eq. (32) for \(Q_{f}=0.01\), \(V_{0}=0.1\). The third curve is pretty close to a rectangular fracture. Further increase of \(b_{2}\) generates shapes with two maxima and a minimum in the middle. Moreover, at even larger \(b_{2}\) the curve ABC (Fig. 2a) intersects the abscissa axis, i.e. the corresponding flow rate function generates physically unrealistic flows, the situation common in the inverse BVP approach (see, e.g. Kacimov [38]).

The complex potential is determined by Eq. (VII.2.13) of [15]:

$$\begin{aligned} w(z)=V{ }_0^{*} z+\frac{1}{2\pi }\int \limits _{-1}^1 {\frac{\omega (u)}{u-z}} \hbox {d}u \end{aligned}$$
(33)

that for the selected truncation of the series (30) gives

$$\begin{aligned} w(z)=(V{ }_0-2Q_f /\pi )^{*} z+Q_f /\pi \left( {1-z^{2}} \right) \log \frac{z-1}{z+1}-\frac{2Q{ }_{f} b_{2}}{3\pi }\int \limits _{-1}^{1} {\frac{(1-u^{2})^{3/2}}{u-z}} \hbox {d}u. \end{aligned}$$
(34)

In Eq. (34), we evaluated the integral by the NIntegrate routine of Mathematica, and used the routines Im and ContourPlot. The obtained streamlines are shown in Fig. 5 for Pilatovskii’s fracture \(Q_{f}=0.01\), \(V_{0}=0.1\), \(b_{2}=0\); curves 1–7 indicate the streamlines \(\psi = 0.0025\), 0.005, 0.01, 0.0125, 0.015, 0.02, 0.03 (curves 1–7, correspondingly).

4 Concluding remarks

We emphasize that the analytical treatment of Robin’s boundary value problems to the Laplace’s 2-D equation is well documented in textbooks on mathematical physics (see, e.g. [27]) and in specialized compendia on, for example, heat conduction (Carslaw and Jaeger [39]) or hydrogeology (Bruggeman [40]). The novelty of our work is in applications to engineering problems of 2-D Darcian flows emerging in hydrofracking, permeable reactive barriers and natural groundwater flows in vicinity of lenticular porous bodies, the permeability of which is much higher than that of the ambient soil or rock. The fractures, which we considered, are “passive”, i.e. there are no sinks (production-injection wells) in the modelled flow domains. Pilatovskii [15] indicated the way how the fractures (not necessarily planar) placed in the formations with “singularities” (wells) can be modelled. Correspondingly, the case of a rectangular (Sect. 2) and arbitrary curvilinear (Sect. 3) fractures placed in ambient fields with sinks-sources can be studied similarly to what we did in this paper.

Optimal shapes of fractures can be sought in the following manner: in the flow rate function, a finite set of coefficients, say, (\(b_{2}, b_{4}\)) can be selected as a dyad of control variables; then-by analogy with aerodynamics ([1])—a fracture of best hydraulic properties can be found. For the Prandlt wings, the lift–drag forces and force moments were criteria in optimization ([1]). For example, for a thin parabolic airfoil it is easy to show that an arc of a fixed length has a maximum lift at a certain degree of bending. In subsurface mechanics, other hydraulic properties are of interest, e.g. travel times of marked particles along streamlines and drops of the piezometric head along them, magnitudes of Darcian velocity in evaluation of hydrodynamics dispersion, etc. (Kacimov and Yakimov [41]). These quantities can be extremized by variation of (\(b_{2}, b_{4}\), \(b_{6}, {\ldots }\)) similarly to what we did for elliptical lenses (Kacimov et al. [17], Kacimov and Obnosov [8]).