2 Analytical solution to a mixed BVP, Poisson”s PDE

In a cross-section perpendicular to the macroscopically vertical direction Of (ascending evaporation), we assume that the pore makes a triangle $OO_{1}O_{2}$ having arbitrary sides $c_{1}$ , $c_{2}$ , $c_{3}$ , (Figure 1(a)). As a geometrical input of our model, we consider the vicinity of one corner (O) of this triangle where water is contained in a “bridge” bounded by two menisci, $AMD$ and $CNB$ , and two pore walls making an angle $\alpha $ , $0 <\alpha <\pi $ (Figure 1(b)). The menisci are circular arcs of constant radii, $r_{1}$ and $r_{2}$ , respectively. We assume that the contact angle $\theta =\pi /2$ .

The gradient of soil water pressure $\nabla p_{w} $ is assumed to be constant, as in [Reference Kacimov, Maklakov, Kayumov and Al-Futaisi22]. Such gradients emerge due to, for example, an increase in the disjoining and capillary pressure caused by the above-mentioned gradual “evaporative thinning” (that is, the decrease of $r_{2}$ and increase of $r_{1}$ in the Of direction) of the water bridges. A higher evaporation rate from menisci of larger diameters (just above the capillary fringe in desert solonchaks), as compared with the vicinity of the super hot desert surface, is well documented.

We introduce a complex physical coordinate $z=x+{i}y=r{e}^{{i}\varphi } $ . The water microfluidic velocity obeys the Poisson equation [Reference Avkhadiev and Kacimov6]: that is

(2.1) $$ \begin{align} \Delta u(r,\varphi )=\frac{\partial ^{2} u}{\partial r^{2} } +\frac{1}{r} \frac{\partial u}{\partial r} +\frac{1}{r^{2} } \frac{\partial ^{2} u}{\partial \varphi ^{2} } =-e,\quad e=\nabla p_{w} /\mu _{w}, \end{align} $$

where $\mu _{w}$ is the viscosity of the pore water. The right-hand side, e, in equation (2.1) is a given constant as a model hydrodynamic input (see [Reference Kacimov, Maklakov, Kayumov and Al-Futaisi22, Reference Philip25] for more details of similar flows in a circular digon and triangle).

The domain $G_{z}$ in Figure 1(c) is a curvilinear tetragon made of two straight segments $DC$ and $BA$ and two arcs $AMD$ and $CND$ . At the boundary of $G_{z}$ , the following conditions hold:

(2.2) $$ \begin{align} \begin{array}{ll} u(r{e}^{\pm {i}\alpha /2} )=0 & r_{1} < r< r_{2} , \\ \dfrac{\partial u(r_{k} {e}^{{i}\varphi } )}{\partial r} =0 & k=1,2,\quad -\alpha /2<\varphi <\alpha /2, \end{array} \end{align} $$

where r is a radial coordinate from the corner O and $\varphi $ is an angular coordinate counted from the axis of symmetry of the bridge, that is, $OMN$ . The pair $r_{1}$ and $r_{2}$ provide another geometrical input of our model. The no-slip (the first line in equation (2.1)) and no-shear-stress (second line) conditions physically mean that the rigid skeleton of our porous medium is static (always true) and $\mu _{a}$ is almost zero, as was assumed in [Reference Kacimov, Maklakov, Kayumov and Al-Futaisi22]. The latter assumption is justified by the two-order difference in viscosities of air and pore water. The main output of the model given by (2.1)–(2.2) is the velocity distribution $u(r, \varphi )$ and ensuing areal integral of this quantity over $G_{z}$ , which gives the volumetric flow rate, and hence the water phase permeability. Extension of the model given by (2.1)–(2.2) to more complex menisci can be found in [Reference Finn15] among other applications to microgravity and soil physics.

The solution of a BVP with the $u = 0$ condition along the whole wall of our $G_z$ was obtained by Sparrow et al. [Reference Sparrow, Chen and Jonsson32], and it was reported by Shah and London [Reference Shah and London31]. Unfortunately, no details of the solution are presented in [Reference Shah and London31, Reference Sparrow, Chen and Jonsson32]. Moreover, equation (394) in [Reference Shah and London31] does not match equation (7) in [Reference Sparrow, Chen and Jonsson32]. Sparrow et al. used the method of separation of variables for solving their BVP for a harmonic function obeying Dirichlet”s boundary conditions, that is, they did not use the theory of holomorphic functions.

Aver $'$ yanov”s [Reference Aver’yanov4] solution to PDE (2.1) can be readily reformulated for $G_{z}$ in the following manner: the boundary conditions are no slip along $DA$ in Figure 1(c) and no shear along $DC$ , $CB$ and $CA$ .

Figure 2 (a) Rectangle $G_{z_1}$ . (b) Reference half-plane $G_{\zeta }$ .

Similar to the work in the papers [Reference Kacimov, Yakimov and Šimůnek21, Reference Kacimov, Maklakov, Kayumov and Al-Futaisi22, Reference Sparrow, Chen and Jonsson32], we reduce the BVP (2.1)–(2.2) for the Poisson equation to a BVP for the Laplace equation. For this purpose, we use the transformation

(2.3) $$ \begin{align} \begin{array}{ll} {z_{1} =\log (z/r_{1} )=x_{1} +{i}y_{1} =\log (r/r_{1} )+{i}\varphi , \quad x_{1} =\log (r/r_{1} ),\quad y_{1} =\varphi ,} \\[.2cm] {\Phi (z)=u(r,\varphi )+er^{2} /4, \quad \Phi _{1} (z_{1} )=\Phi (r_{1} \exp \, z_{1} ),} \end{array} \end{align} $$

which makes the function u in equation (2.1) harmonic in a rectangle $G_{z_{1} } $ (Figure 2(a)): that is,

(2.4) $$ \begin{align} \Delta \Phi _{1} (x_{1} ,y_{1} )=0 \quad \textrm{for} \quad z_{1} =x_{1} +{i}y_{1} \in G_{z_{1}}. \end{align} $$

The boundary conditions (2.2), (2.3) in $G_{z_{1} } $ become

(2.5) $$ \begin{align} AB, DC: &\, \Phi _{1} (x_{1} ,\pm \alpha )=0.25\mathrm{\; }e\mathrm{\; }r_{1}^{2} {e}^{2x_{1} } \quad 0<x_{1} <l=\log (r_{2} /r_{1} ), \nonumber\\ AD: & \,\dfrac{\partial \Phi _{1} (0,y_{1} )}{\partial x_{1} } =\frac{e\mathrm{\; }r_{1}^{2} }{2} \quad -\alpha /2<y_{1} <\alpha /2, \\ BC:& \, \dfrac{\partial \Phi _{1} (l,y_{1} )}{\partial x_{1} } =\frac{e\mathrm{\; }r_{2}^{2} }{2} \quad -\alpha /2<y_{1} <\alpha /2.\nonumber \end{align} $$

Clearly, the function $\Phi _{1} (x_{1} ,-y_{1} )$ is a solution of the mixed BVP (2.4) and (2.5) together with $\Phi _{1} (x_{1} ,y_{1} )$ . Hence,

$$ \begin{align*}\Phi _{1} (x_{1} ,-y_{1} )=\Phi _{1} (x_{1} ,y_{1} ),\end{align*} $$

due to the uniqueness theorem [Reference Henrici18].

Next, we introduce the complex potential function

(2.6) $$ \begin{align} \Omega _{1} (z_{1} )=\Phi _{1} (x_{1} ,y_{1} )+{i}\Psi _{1} (x_{1} ,y_{1} ), \end{align} $$

which is holomorphic in $G_{z_{1} } $ . Here the function $\Psi _{1} (x_{1} ,y_{1} )$ , being complex conjugate to $\Phi _{1} (x_{1} ,y_{1} )$ , satisfies the symmetry condition

$$ \begin{align*} \Psi _{1} (x_{1} ,-y_{1} )=-\Psi _{1} (x_{1} ,y_{1} ). \end{align*} $$

Due to equations (2.5) and the Riemann–Schwartz symmetry principle [Reference Gakhov17], the function (2.6) admits analytic continuations through the vertical sides of the boundary of the rectangle $G_{z_{1} } $ . Hence, the Cauchy–Riemann conditions (CRC) hold up to these sides. Using the CRC $\partial \Phi _{1} /\partial x_{1} =\partial \Psi _{1} /\partial y_{1} $ and the last two conditions of (2.5) and integrating, we get

(2.7) $$ \begin{align} \Psi _{1} (0,y_{1} )=0.5er_{1}^{2} y_{1} ,\quad \Psi _{1} (l,y_{1} )=0.5er_{2}^{2} y_{1} \quad -\alpha <y_{1} <\alpha. \end{align} $$

It is noteworthy that an arbitrary constant arising after integrating the CRC condition is equal to zero due to the above-mentioned symmetry property of the function $\Psi _{1} (x_{1} ,y_{1} )$ .

We map conformally $G_{z_{1} } $ onto the upper half-plane $\mathrm {Im}\zeta>0$ , $G_{\zeta }={\mathbb C}^+ $ (Figure 2(b)), with the correspondence of the points

$$ \begin{align*}C\to -1/\lambda , \quad D\to -1, \quad A\to 1, \quad B\to 1/\lambda , \end{align*} $$

where $0<\lambda <1$ . The Schwarz–Christoffel formula [Reference Polubarinova-Kochina26] gives

(2.8) $$ \begin{align} z_{1} (\zeta )=-\frac{{i}\alpha }{2K(\lambda )} \int _{0}^{\zeta }\frac{{d}t}{\sqrt{(1-t^{2} )(1-\lambda ^{2} t^{2} )} } =-\frac{{i}\alpha }{2K(\lambda )} \mathrm{F}(\arcsin \zeta ,\lambda ) , \end{align} $$

where $\mathrm {F}(\arcsin \zeta ,\lambda )$ and $K(\lambda )$ are, respectively, incomplete and complete elliptic integrals of the first kind (see [Reference Abramowitz and Stegun1, formulae 17.2.7, 17.3.1]). The modulus $\lambda $ is determined from the condition

(2.9) $$ \begin{align} K'(\lambda )/K(\lambda )=2l/\alpha, \end{align} $$

where $K'(\lambda )=K(\lambda ')$ , $\lambda '=\sqrt {1-\lambda ^{2} } $ . We use the FindRoot routine of Mathematica [Reference Wolfram36] and find $\lambda $ from equation (2.9).

Combining conditions (2.5) and (2.7), we obtain the following mixed BVP for the function $\Omega (\zeta )=\Omega _{1} (z(\zeta ))$ in $G_{\zeta } $

(2.10) $$ \begin{align} \begin{array}{ll} \mathrm{Re}\Omega (\xi )=0.25er_{1}^{2} {e}^{2x_{1} (\xi )} & 1<|\xi |<1/\lambda , \\[.2cm] \mathrm{Im}\Omega (\xi )=0.5er_{1}^{2} y_{1} (\xi ) & -1<\xi <1, \\[.2cm] \mathrm{Im}\Omega (\xi )=0.5er_{2}^{2} y_{1} (\xi ) & |\xi |>1/\lambda. \end{array} \end{align} $$

Here $x_{1} (\xi )$ and $y_{1} (\xi )$ are real and imaginary parts of the function (2.8), respectively: that is,

(2.11) $$ \begin{align} \begin{aligned} x_{1} (\xi )&=\dfrac{\alpha }{2K(\lambda )} \int_{1}^{\xi }\frac{{d}t}{\sqrt{(t^{2} -1)(1-\lambda ^{2} t^{2} )} } =\frac{\alpha }{2K(\lambda )} F\bigg(\!\arcsin \dfrac{\sqrt{1-\xi ^{-2} } }{\lambda '} ,\lambda '\bigg) \quad 1<\xi <1/\lambda , \\ y_{1} (\xi )&=-\dfrac{\alpha }{2K(\lambda )} \int_{0}^{\xi }\dfrac{{d}t}{\sqrt{(1-t^{2} )(1-\lambda ^{2} t^{2} )} } =\dfrac{\alpha }{2K(\lambda )} F(\arcsin \xi ,\lambda ) \quad-1<\xi <1, \\ y_{1} (\xi )&=\dfrac{\alpha }{2K(\lambda )}\! \int_{1/(\lambda \xi )}^{1}\!\dfrac{{d}t}{\sqrt{(1-t^{2} )(1-\lambda ^{2} t^{2} )} } -\dfrac{\alpha }{2} \\ &=\frac{\alpha }{2K(\lambda )} F\bigg(\!\arcsin \sqrt{\dfrac{\xi ^{2} -1/\lambda ^{2} }{\xi ^{2} -1} } ,\lambda \bigg)-\dfrac{\alpha }{2} \quad \xi>1/\lambda . \end{aligned} \end{align} $$

It is clear that $x_{1} (\xi )=x_{1} (-\xi )$ for $-1/\lambda <\xi <-1$ and $y_{1} (\xi )=-y_{1} (-\xi )$ for $\xi <-1/\lambda $ .

Thus, we have to find $\Omega (\zeta )$ , which is holomorphic in $G_{\zeta } $ . This function must satisfy boundary conditions (2.10). It must be bounded at all transition points $\pm 1$ , $\pm 1/\lambda $ and at infinity in Figure 2(b). The index of the BVP (2.10) in the indicated class of functions is $-2$ , and, generally speaking, one solvability condition must be satisfied [Reference Gakhov17]. As will be shown below, this solvability condition takes place automatically due to the properties of functions (2.11).

Solution of problem (2.10) bounded at transition points, but generally, not at infinity, [Reference Henrici18] is

$$ \begin{align*} \Omega (\zeta )=\frac{er_{1}^{2} \omega \, (\zeta )}{2\pi } \bigg(\!\int_{-1}^{1}\frac{y_{1} (\xi )}{\omega \, (\xi )} \frac{{d}\xi }{\xi -\zeta } +\int_{L_\lambda }\frac{{e}^{2x_{1} (\xi )} }{2|\omega \, (\xi )|} \frac{{d}\xi }{\xi -\zeta } +{e}^{2l} \int_{1/\lambda }^{-1/\lambda }\frac{y_{1} (\xi )}{\omega \, (\xi )} \frac{{d}\xi }{\xi -\zeta } \bigg),\end{align*} $$

where $L_\lambda =(-1,-1/\lambda )\cup (1,1/\lambda )$ , and the last integral is taken from $1/\lambda $ to $-1/\lambda $ through infinity. The branch of the function

(2.12) $$ \begin{align}\omega (\zeta )=\sqrt{(1-\zeta ^{2} )(1-\lambda ^{2} \zeta ^{2} )} , \end{align} $$

is fixed in the upper half-plane by the condition of its positivity in the interval (-1,1). This branch is real and negative at $\zeta =\xi \in (-\infty ,-1/\lambda )\cup (1/\lambda ,\infty )$ . The branch is pure imaginary with negative and positive imaginary parts at the intervals $(1,1/\lambda )$ and $(-1/\lambda ,-1)$ , respectively. It may seem that the $\Omega (\zeta )$ obtained in equation (2.12) is unbounded at infinity. If, however, we take into account the properties of functions (2.11) and (2.12), this solution can be rewritten as

(2.13) $$ \begin{align} \Omega (\zeta )=\frac{er_{1}^{2} \omega (\zeta )}{\pi } \bigg(\!\int_{0}^{1}\frac{\xi y_{1} (\xi )}{\omega (\xi )} \frac{{d}\xi }{\xi ^{2} -\zeta ^{2} } +\frac{1}{2} \int_{1}^{1/\lambda }\frac{\xi {e}^{2x_{1} (\xi )} }{|\omega (\xi )|} \frac{{d}\xi }{\xi ^{2} -\zeta ^{2} } +{e}^{2l} \int_{1/\lambda }^{\infty }\frac{\xi y_{1} (\xi )}{\omega (\xi )} \frac{{d}\xi }{\xi ^{2} -\zeta ^{2} } \bigg). \end{align} $$

From equation (2.13), it is clear that our solution is bounded at $\zeta \to \infty $ , as it should be. Note that the functions (2.8), (2.12) and (2.13) satisfy the same symmetry condition

(2.14) $$ \begin{align}z_{1} (\zeta )\equiv \overline{z_{1} (-\bar{\zeta })}, \quad \omega (\zeta )\equiv \overline{\omega (-\bar{\zeta })}, \quad \Omega (\zeta )\equiv \overline{\Omega (-\bar{\zeta })}. \end{align} $$

The total flow rate through $G_{z}$ , due to equation (2.3), is

$$ \begin{align*} Q=\int_{G_{z} }u(x,y)\,{d}s =\int _{G_{z} }\Phi (x,y)\,{d}s -\frac{e\alpha }{16} (r_{2}^{4} -r_{1}^{4} ). \end{align*} $$

The last integral can be transformed as

$$ \begin{align*} Q_{0} =\int _{G_{z} }\Phi (x,y)\,{d}s =\frac{{i}}{2} \int _{G_{z} }\Phi (x,y)\,{d}z\wedge {\kern 1pt} {d}\bar{z}=\frac{{i}r_{1}^{2} }{2} \int _{G_{z_{1} } }\Phi _{1} (x_{1} ,y_{1} )e^{2x_{1} } \,{d}z_{1} \wedge {\kern 1pt} {d}\bar{z}_{1} , \end{align*} $$

where $\wedge $ is the symbol for the exterior or wedge product. Using the transformation (2.8), we obtain

(2.15) $$ \begin{align} Q_{0} =\frac{{i}r_{1}^{2} \alpha ^{2} }{16K^{2} } \int_{{\mathbb C}^{+} }\frac{\Omega (\zeta )+\overline{\Omega (\zeta )}}{\omega (\zeta )\overline{\omega (\zeta )}} e^{2\mathrm{Re} z_{1} (\zeta )} \,{d}\zeta \wedge {d}\bar{\zeta }=\frac{r_{1}^{2} \alpha ^{2} }{4K^{2} } \int_{-\infty }^{\infty }\int_{0}^{\infty }\frac{\mathrm{Re}\Omega (\xi +{i}\eta )}{|\omega (\xi +{i}\eta )|^{2} } e^{2x_{1} (\xi +{i}\eta )} \,{d}\eta\, {d}\xi. \end{align} $$

Using the Stokes formula, $Q_{0}$ is transformed into a contour integral along the real axis

$$ \begin{align*}Q_{0} =\frac{{i}r_{1}^{2} \alpha ^{2} }{16K^{2} } \bigg(\!\int_{-\infty }^{\infty }\frac{e^{\overline{z_{1} (\xi )}} }{\overline{\omega (\xi )}} \int_{0}^{\xi }\frac{\Omega (t)e^{z_{1} (t)} }{\omega (t)} {d}t\,{d}\xi -\int_{-\infty }^{\infty }\frac{e^{z_{1} (\xi )} }{\omega (\xi )} \int_{0}^{\xi }\frac{\overline{\Omega (t)}e^{\overline{z_{1} (t)}} }{\overline{\omega (t)}} {d}t\,{d}\xi \bigg).\end{align*} $$

Finally, using the symmetry properties (2.14), we obtain

(2.16) $$ \begin{align} Q=-\frac{r_{1}^{2} \alpha ^{2} }{4K^{2} } \mathrm{Im}\int_{0}^{\infty }\frac{e^{\overline{z_{1} (\xi )}} }{\overline{\omega (\xi )}} \int_{0}^{\xi }\frac{\Omega (t)e^{z_{1} (t)} }{\omega (t)} {d}t\,{d}\xi -\frac{e\alpha }{16} (r_{2}^{4} -r_{1}^{4} ). \end{align} $$

We used the NIntegrate routine of Wolfram”s Mathematica [Reference Wolfram36] to compute the integrals in equations (2.15) and (2.16). It turned out that computation Q using equation (2.16) required up to three days (for example, on a PC Intel(R) Core(TM) i7-4770 CPU 3.40GHz, 3401 IHz, cores: 4, logical processors: 8, which we used) for each value of $r_{1} $ . If one uses equation (2.15) (as we did) the computations take up to an hour.

We introduce the dimensional quantities $Q^{*}=Q/( e r_{2}^{4})$ and $r^{*}=r_{1}/r_{2}$ . Table 1 presents the values of $Q^{*}$ for $\alpha =\pi /3$ and five values of $r^{*}$ .

Table 1 Dimensionless flow rate $Q^{*}$ computed using equations (2.15) and (2.16) for tetragons $G_{z}$ having $\alpha =\pi /3$ and five values of $r^{*}$ .

Remark 1 In addition, for simplification of the expression $\Omega (\zeta )$ in equation (2.15), we used the substitution $\xi \to -1/(\lambda \xi )$ and transformed the last term in equation (2.13) to give

$$ \begin{align*}{e}^{2l} \int_{1/\lambda }^{\infty }\frac{\xi y_{1} (\xi )}{\omega (\xi )} \frac{{d}\xi }{\xi ^{2} -\zeta ^{2} } ={e}^{2l} \int_{-\infty }^{-1/\lambda }\frac{\xi y_{1} (\xi )}{\omega (\xi )} \frac{{d}\xi }{\xi ^{2} -\zeta ^{2} } =-\lambda {e}^{2l} \int_{0}^{1}\frac{\xi y_{1} (\xi )}{\omega (\xi )} \frac{{d}\xi }{1-\lambda ^{2} \xi ^{2} \zeta ^{2} } .\end{align*} $$

To obtain this result, we used the identities

$$ \begin{align*}\omega (-1/(\lambda \zeta ))\equiv -\omega (\zeta )/(\lambda \zeta ^{2} ),\quad z_{1} (-1/(\lambda \zeta ))\equiv l-z_{1} (\zeta ). \end{align*} $$

3 Other applications of water flow at high $S_{w}$ in the pore core with static air in the corners

A curvilinear tetragon, as considered in Section 2 for the case of “thin bridges” and small $S_{w}$ , can be used as a generic mobile phase conduit in the opposite limit of large $S_{w}$ (close to 1). Indeed, let viscous water occupy almost the whole pore (see an equilateral triangle in Figure 3(a)) and move in the Of direction (not shown in Figure 3). Tiny pockets of immobile air are “entrapped” near the three corners in Figure 3(a) (see [Reference Blunt9]) which, in applications to reservoir engineering, can be considered (at the REV-scale) as an early (primary) stage of crude oil recovery with a minor proportion of gas formation in the rock. The water flow domain in Figure 3(a) is a circular sextagon, with point E being its centre. Owing to symmetry, we can consider the curvilinear pentagon $ADE_{C}EE_{B}A$ , that is, the third part of $G_{z}$ in Figure 3(a). Along the broken line $E_{C}EE_{B}$ the shear stress is zero (owing to symmetry). We can approximate this curvilinear pentagon $DE_{C}EE_{B}A$ by a curvilinear tetragon $DE_{C}NE_{B}A$ , in which the wedge $E_{C}EE_{B}$ is replaced by a circular arc $E_{C}NE_{B}$ , that is, we arrive at $G_{z}$ in our Figure 1. Consequently, in $ADE_{C}NE_{B}A$ of Figure 3, all analysis in Section 2 is valid.

Figure 3 (a) Flow of water at high $S_{w}$ in a cylindrical pore channel, the cross-section of which is an equilateral triangle with immobile air pockets entrapped near the corners. (b) Flow of air at low $S_{w}$ in a cylindrical pore channel, the cross-section of which is an equilateral triangle with immobile water pockets entrapped near the corners.

In Figure 3(b), another limiting case of very low $S_{w}$ is depicted. Now the pore gas is moving in the Of direction. The right-hand side in equation (2.1) becomes $e=\nabla p_{a} /\mu _{a} $ , where $\mu _{a}$ is the viscosity of the air and $\nabla p_{a} $ is the gradient of the gas pressure along the pore. Water in Figure 3(b) is entrapped near the triangle corners, that is, it is immobile. Similarly to the case in Figure 3(a), we make a circular tetragon with mixed boundary conditions on its sides such that the analysis in Section 2 is valid. Applications of this pore-scale model are to the emission of in-pore CO ${}_{2}$ from the root zone, which has much higher concentrations of CO ${}_{2}$ than the superjacent atmospheric air. Evaluation of the fluxes of gases through the vadose zone are of utmost importance (see [Reference Ben-Noah and Friedman8, Reference Scanlon, Nicot, Massmann and Warrick30]), especially in the context of the impact of allegedly deleterious greenhouse gases (generated by rhizomes in the soil) on the human race and terrestrial ecosystems.

4 Concluding remarks on limitations and perspectives of our model

Our model and machinery of solving mixed BVPs for holomorphic functions allows us to use a combination of slip and no-slip boundary conditions, which is applicable when one phase is an immobile gas while the pore is polygonal rather than the Aver”yanov variably wet circle. In other words, we extend the cases of digon and trigon with mixed boundary conditions, considered in [Reference Kacimov, Maklakov, Kayumov and Al-Futaisi22, Reference Philip25], to the case of a circular tetragon, on the opposite sides of which these conditions are imposed. The motion of moisture in a “bridge”, which is modelled by such a tetragon, is described by Poisson”s equation. With the help of the Keldysh–Sedov-type formula [Reference Henrici18], we obtained an integral representation for the in-pore velocity.

We have doubts about the solution obtained by Sparrow et al. [Reference Sparrow, Chen and Jonsson32] to the BVP with $u= 0$ on the whole boundary of the tetragon in Figure 1(c). The lack of comparisons in [Reference Sparrow, Chen and Jonsson32] (as well as in [Reference Shah and London31]) with the earlier explicit solution to the Poisson equation in a circular triangle, obtained by Saint–Venant [Reference Timoshenko and Goodier34] is confusing. Indeed, the solution obtained in [Reference Sparrow, Chen and Jonsson32] should degenerate into Saint–Venant”s solution when the small-radius meniscus in Figure 1 vanishes, which was not demonstrated by Sparrow et al. [Reference Sparrow, Chen and Jonsson32].

The estimates of the flow rate Q through polygonal cross-sections with no-shear menisci can be attempted with the help of isoperimetric inequalities [Reference Pólya and Szegö27].

Recent video recording techniques of microfluidics processes can be used as experimental tools for measuring the in-pore parameters $\alpha $ , $r_1$ , $r_2$ in Figure 1.

Solute transport on the microscale was not studied in our paper. In future, the advective dispersion equation on the macroscale (see [Reference Aver’yanov5]) can be downscaled to the pore level. The root zone salinization, caused by evapotranspiration from shallow water tables, is a daunting problem in Australia and other arid countries, where applied mathematicians may, like Aver $'$ yanov, transcend from microfluidics to the scale of engineered drainage of large salinized crop fields.

The limitations of our flow model in Section 2 are as follows.

1. (1) Viscosity of water is constant (in reality, it varies along the Of-axis, because the soil temperature of deserts has a strong diurnal gradient in the vertical direction).

2. (2) The flow domain $G_{z }$ does not vary in the Of direction (in reality, water evaporates into the pore air and $G_{z }$ dwindles from the capillary fringe to the topsoil).

3. (3) The water pressure gradient is constant (in reality, it is not because the decrease of $r_{1}$ and $r_{2}$ ) in the Of direction causes variation of this gradient.

4. (4) The pores are cylindrical and flow is two-dimensional (in reality, the pores are not cylindrical and flow is microscopically three-dimensional).

5. (5) The contact angle is $90^{\circ }$ (in reality, most soils are either hydrophilic or hydrophobic).

6. (6) In Section 2, air is immobile (in reality, both water and air move, and not always co-currently).

7. (7) Water “bridges” are assumed to be stable to potential perturbations (in reality, the menisci in Figure 1 may break into Avery $'$ anov-type films that simply “coat” the wall of the polygonal pore, which gives the case of flow considered by Kacimov et al. [Reference Kacimov, Maklakov, Kayumov and Al-Futaisi22]).

8. (8) Gravity is ignored, that is, capillarity dominates the fluid behaviour, along with the hydrodynamic pressure gradient along Of (in thick vadose zones, the gravity force is explicitly counted on the right-hand side of the Poisson equation).

9. (9) Menisci are circular arcs (in reality, they are not [Reference Finn15]).