Experimental modelling of hydraulic parameters for fluid flow in stratified porous media

Abstract

From works in in-situ seepage through dams and laboratory experiments using Layered Heterogeneous Porous Media (LHPM), it has been noted that a refraction-like phenomenon, such as that experienced in light, affects fluids’ flowlines when crossing the contact interface of layers characterised by different porosity viz-a-viz permeability. This concept has many applications in fluid dynamics, such as the dispersion process in stratified media. Currently, no study exists that models and analyses the relationship between the porosities of two layers in contact and the resulting flowline refraction using validated LHPM data. Hence, this work aims to establish a relationship between the porosity ratio $\left({\mathrm{\Phi }}_r\right)$ of stratified media made up of two layers and the refraction angle $\left({\theta }_{\mathit{max}}\right)$ of the maximum volume flux $\left(q_{\mathit{max}}\right)$. Volume flux data from the flow of a single-phase fluid through five LHPMs with $\underset{\mathit{r}}{\mathrm{\Phi }}$ ranging from 0.8325 to 0.9524 were used in the modelling. The flow was oriented from the lower to the higher porosity vis-à-vis permeability layer. It was found that ${\theta }_{\mathit{max}}$ which is also the refraction angle of the peak solute plume flux, refracts away from the normal as ${\mathrm{\Phi }}_r$ reduces. This indicates an increase in the dilution rate vis-à-vis spread of plumes with the reduction in homogeneity between the two layers. Also, ${\theta }_{\mathit{max}}$ does not correlate with the stratification inclination $\left(\mathit{\alpha }\right)$, but the $q_{\mathit{max}}$ which is also the peak solute plume flux correlates with $\mathit{\alpha }$. Furthermore, an efficient model, which is the best-unbiased estimator, with $R^2=0.99$ was derived. Findings from this work can help better understand solute plume dispersion and the general fluid flow dynamics in stratified media such as capillary barrier effect covers for pollution control and hydrocarbon reservoirs.

Keywords:

• capillary barrier effect
• layered porous media
• porosity
• preferential flow
• dispersion

1. Introduction

When fluid flow through any medium, rocks and soil inclusive, its flow path and rate are influenced by several petrophysical properties such as tortuosity, pore size distribution, cracks, fissures, porosity, permeability, and inclination, which determine volume flux (q) at every instant, hence making in-situ flow dynamics complex (Wei et al., 2021; Alabi & Sedara, 2020; H. H. Wang et al., 2020; X. Li et al., 2020; Sun & Zhang, 2020; Govindarajan, 2019; Cao et al., 2016). Therefore, Intermediate Scale Experiments (ISEs) carried out in the laboratory are used to study the specific properties of interest, e.g., Preferential Flow Features (PFF) and Barrier Flow Features (BFF) and how they affect the flow dynamics (Solnyshkina et al., 2021; X. Li et al., 2020; Sebben & Werner, 2016b; Silliman et al., 1998). One of such ISEs is Layered Heterogeneous Porous Media (LHPM), a proxy used to study the effect of stratifications on flow dynamics (Alabi & Akanni, 2021, 2022a; Ku et al., 2021; X. Li et al., 2020; Liang et al., 2019). LHPMs have found application in geotechnics, CO₂ sequestration, enhanced recovery of oil and gas, landfills’ leachate collection systems, and pollution control systems (Alabi & Akanni, 2022b; Banihashem & Karrabi, 2020; Hosking et al., 2018; Hu et al., 2020; March et al., 2018; Sai et al., 2020; Salimzadeh & Khalili, 2014; Zhang & Yuan, 2019).

From works in in-situ seepage through dams and laboratory experiments using LHPMs, it has been noted that a refraction-like phenomenon, such as that experienced in light, affects fluids’ flowlines when crossing the contact interface of layers characterised by different porosity viz-a-viz permeability (Casagrande, 1937; Cedergren, 1997; Eagleman & Jamison, 1962). However, unlike Snell’s law which is a sine law, groundwater refraction obeys tangential law (Sebben & Werner, 2016a).

Figure 1b illustrates the refraction schematic diagram for an LHPM of increasing layer hydraulic conductivity (K) in the flow direction. The flow originates from the top layer (L1) of $\mathit{K}=K_1$, into the bottom layer (L2) of $\mathit{K}=K_2$, where $K_1<K_2$. This flow regime is also known as Ascending Flow (AF) system. In such a system, the flowline refracts from the normal in L2 (Alabi & Akanni, 2022a). The observed refraction is due to the increased fluid velocity in L2 as the K increases (Pavlovskaya et al., 2018; Sebben & Werner, 2016b). It results in dispersion towards a higher refraction angle from the normal because the fluid flow field shows the instantaneous and free stream velocity (Ferdows et al., 2022; Sebben & Werner, 2016b). Figure 2 illustrates the reverse flow, i.e. flow from *L1 into *L2, where ${ }^X{K}_1{>}^X{K}_2$. This flow system is also known as Descending Flow (DF). In the DF system, the flowline refracts towards the normal in *L2 (Alabi & Akanni, 2022a). Since this study focuses on the AF system, more focus would be placed on it.

Figure 1. a) Flownet showing groundwater refraction in a heterogeneous medium. Areas with the same K have equal-sized “squares.” b) A simplified flownet refraction in an AF system (heterogeneous isotropic medium). The flow is from L1 into L2. It is characterised by two flowlines (yellow dash lines), two different equipotential lines (red dash lines), ${\mathrm{\Phi }}_1$ as porosity for L1, $K_1$ as hydraulic conductivity for L1, ${\mathrm{\Phi }}_2$ as porosity for L2, and $K_2$ as hydraulic conductivity for L2 (Adapted from Bos, 2006; Freeze & Cherry, 1979; Leliavsky, 1955).

Figure 2. Refraction in a DF system (heterogeneous isotropic medium). The flow is from *L1 into *L2. It is characterised by two flowlines (yellow dash lines), two different equipotential lines (red dash lines), ${ }^X{\mathrm{\Phi }}_1$ as porosity for L1, ${ }^X{K}_1$ as hydraulic conductivity for *L1, ${ }^X{\mathrm{\Phi }}_2$ as porosity for *L2, and ${ }^X{K}_2$ as hydraulic conductivity for L2 (Adapted from Bos, 2006; Freeze & Cherry, 1979; Leliavsky, 1955).

From Figure 1b, the tangential law of groundwater refraction can be derived by first considering the mass continuity equation. It states that the discharge in both layers must be equal irrespective of the K value of each layer (Bos, 2006; Y. Li et al., 2022). Therefore,

(1) $\mathrm{\Delta q}=i_1{K}_1\mathit{a}cos{\theta }_1=i_2{K}_2\mathit{a}cos{\theta }_2$(1)

where $\mathrm{\Delta }\mathit{q}$ is the volumetric discharge that flows between the flowlines (F1 and F2), $i_1$ is the hydraulic gradient of L1, is the hydraulic conductivity of L1, $\mathit{a}$ is the unit area between two fluid particles following two different flowlines (F1 and F2), ${\theta }_1$ is the entry angle, $i_2$ is the hydraulic gradient of L2, $K_2$ is the hydraulic conductivity of L2, and ${\theta }_2$ is the refraction angle.

The i for the two layers can be expressed as (Bos, 2006),

(2) $i_1=\frac{E_1-E_2}{\mathit{a}sin{\theta }_1};i_2=\frac{E_1-E_2}{\mathit{a}sin{\theta }_2}$(2)

where $E_1$ is equipotential line 1, $E_2$ is equipotential line 2, and all other notations retain their meaning.

Substituting Eqn. 2(2) $i_1=\frac{E_1-E_2}{\mathit{a}sin{\theta }_1};i_2=\frac{E_1-E_2}{\mathit{a}sin{\theta }_2}$(2) into Eqn. 1(1) $\mathrm{\Delta q}=i_1{K}_1\mathit{a}cos{\theta }_1=i_2{K}_2\mathit{a}cos{\theta }_2$(1) results in

(3) $K_1\frac{E_1-E_2}{\mathit{a}sin{\theta }_1}\mathit{a}cos{\theta }_1=K_2\frac{E_1-E_2}{\mathit{a}sin{\theta }_2}\mathit{a}cos{\theta }_2$(3)

Solving Eqn. 3(3) $K_1\frac{E_1-E_2}{\mathit{a}sin{\theta }_1}\mathit{a}cos{\theta }_1=K_2\frac{E_1-E_2}{\mathit{a}sin{\theta }_2}\mathit{a}cos{\theta }_2$(3) , using trigonometry results in the tangential law for groundwater flow (Bos, 2006; Lu et al., 2013);

(4) $\frac{K_1}{K_2}=\frac{tan{\theta }_1}{tan{\theta }_2}$(4)

where all notations retain their meaning.

Furthermore, recent studies have explored the flow dynamics associated with the contact interface in stratifications. Dejam and Hassanzadeh (2022), using analytical and numerical methods, reported the presence of tensorial dispersion and advection at the contact interfaces of LHPMs characterised by different Darcy numbers $\left(D_a\right)$. They reported that such dispersion tensor and advection become symmetric once hydraulic parameters like porosity, thickness, and molecular diffusion become equal, i.e. when the LHPM becomes isotropic and homogeneous. Using numerical analysis, Sebben and Werner (2016b) found that as the L2 width increases, dispersion also increases in AF systems. Kurasawa et al. (2020) reported that for a similarly scaled heterogeneity in LHPMs and randomly porous media, the longitudinal dispersivity of the former would always be more significant. This dispersivity could be enhanced nonmonotonically with the aid of velocity fluctuations as a stochastic spike phenomenon, where the statistically obtainable maximum dispersion is at a correlation length ${\lambda }_{\mathit{max}}\approx 0.25\mathit{y}$, where y is the total layer thickness (Xu & Tartakovsky, 2017). Chen et al. (2022), from ISEs and numerical models, reported the presence of directionally dependent asymmetric Breakthrough Curves (BTCs) in LHPMs with sharp and transition interfaces. The transition interfaces were reported to be associated with more pronounced BTCs, due to more accumulation at such interfaces. These transition interfaces are known to occur naturally or from man-made activities like contamination from nanoparticle-enhanced drilling fluids (L. Bég et al., 2018; L. Wang et al., 2022). Chen et al. (2022) further reported the Mobile and Immobile (MIM) model as being optimally better than the Advection-Dispersion Equation (ADE) model in capturing the directionally dependent transport of solutes. Y. Li et al. (2022) extended the works of (Sebben & Werner, 2016a) and (Sebben & Werner, 2016b) to non-Darcian dispersion of solute plumes in the PFF (faults and fractures) using the finite-element model COMSOL Multiphysics. Such dispersion occurs due to variable porosity and turbulent flow (Usman, Shaheen et al., 2022). Y. Li et al. (2022) reported a 160% increase in solute concentration compared to Darcian models due to the non-Darcian dispersion, resulting from a reduction in interface refraction and specific discharge. This reduced refraction implies that Eqn. 4(4) $\frac{K_1}{K_2}=\frac{tan{\theta }_1}{tan{\theta }_2}$(4) is invalid in non-Darcian conditions. Y. Li et al. (2022) also found that non-Darcian effects are generally negligible in systems where PFFs are filled with porous media, except in rare situations.

All the earlier reviewed works of literature used numerical models to visualise the interface’s refraction and associated solute plume dispersion. For instance, Chen et al. (2022) made use of the Galerkin finite-element software package COMSOL Multiphysics (COMSOL Inc., Burlington, MA, USA), while Dejam and Hassanzadeh (2022) made use of GBAND solver for the resultant fully implicit Finite Difference (FD) discretised equation. Limited literature has reported observations from ISEs, designed and constructed to prove Eqn. 4(4) $\frac{K_1}{K_2}=\frac{tan{\theta }_1}{tan{\theta }_2}$(4) (the tangential groundwater law of refraction) based on the experimental results alone. Of the limited literature, we intend to focus on Popoola et al. 2010a, Popoola et al. 2010b and Popoola et al. 2009, whose full ISE results we have been granted access. The authors reported that the refraction angle $\left(\mathit{\theta }\right)$ of the flowline is associated with the peak volume flux ($q_{\mathit{max}})$, denoted by ${\theta }_{\mathit{max}}$ deflects from the normal for a single-phase fluid flowing from a lower to higher porous vis-à-vis permeable layer (AF system). The reported result aligns with Sebben and Werner (2016a), and it validates the ISE setup as being in line with Eqn. 4(4) $\frac{K_1}{K_2}=\frac{tan{\theta }_1}{tan{\theta }_2}$(4) . Figure 1b illustrates the reported flow dynamics, where ${\theta }_2={\theta }_{\mathit{max}}$.

More recent works have also been done using the ISE’s data (data obtained from Popoola et al., 2010a, Popoola et al., 2010b and Popoola et al., 2009). These subsequent works have further validated the ISE setup and its data. The works show that other flow dynamics deductions from the ISE’s data align with the findings of different numerical studies on AF systems. Using part of the data from the ISE setup as input parameters for 3D numerical analysis, Alabi and Akanni (2021) reported AF systems’ interfaces as points of flow convergence. The reported result agrees with other authors like Ng et al. (2015) and Aubertin et al. (2009), whose studies were based on different data sources for AF system. Alabi and Akanni (2021) further reported the suitability of the LHPM arrangement as capillary barrier effect (CBE) covers for groundwater protection, which aligns with the result of Huo et al. (2022). Furthermore, Alabi & Akanni, (2022a) inferred from applying statistical analysis on the ISE’s data and DF system data (Alabi, 2015; Alabi & Olaleye, 2016) that irrespective of the flow direction, the higher the difference between the hydraulic properties of the layers in contact, the higher the flow turbulence induced on a flowline crossing the interface. This aligns with the numerical simulation result Cherblanc et al. (2007) reported.

Generally, flowline refraction is comparable to solute plume dispersion only in cases where mechanical dispersion and solute molecular diffusion can be ignored, such as cases of single-phase flow being studied (Sebben & Werner, 2016b). Considering that the fluid used in the experiment is single-phase, the ${\theta }_{\mathit{max}}$ can also be referred to as the refraction angle of the peak solute plume flux. This is because, since the fluid is of single-phase type, the $q_{\mathit{max}}$ gives the most volume of the solute plume passing through a unit area of the interface per unit of time. Hence, its refraction $\left({\theta }_{\mathit{max}}\right)$ dictates the dispersion magnitude (Dawe & Marcelle-de Silva, 2005; Govindarajan, 2019). Understanding such dispersion dynamics would offer vital benefits such as a better understanding of nanofluids flooding processes for increasing Enhanced Oil Recovery (EOR) efficiency and dilution processes in other stratified media (Dejam & Hassanzadeh, 2022; Kurasawa et al., 2020; Magnusson et al., 2004; Shamshuddin et al., 2019).

Using the reviewed literature, the LHPM data for this study have been validated. Also, the data have been extended to include dispersion dynamics. To our knowledge, the available works have presented findings on the volume flux and solute plume dispersion dynamics in stratified media using numerical methods based on hydraulic properties like Darcy number and K of the layers in contact (Dejam & Hassanzadeh, 2022; Sebben & Werner, 2016a, 2016b). However, findings based on the porosities of the layers in contact have not been reported from purely experimental data. In addition, no model has been reported that relates the porosities to ${\theta }_{\mathit{max}}$ as a measurement parameter for volume flux and solute plume dispersion dynamics. Therefore, in this work, an analysis and model between the porosities of the layers and ${\theta }_{\mathit{max}}$ would be presented. This could have applications in many areas, such as increasing the efficiency of the flooding processes in EOR, optimising CBE covers, and improving the knowledge of solute dispersal dynamics in other stratified media.

We aim to analyse and model the relationship between the ${\mathrm{\Phi }}_r$ of a two-layer stratified medium and the ${\theta }_{\mathit{max}}$. In addition, we extend the ${\theta }_{\mathit{max}}$ to include dispersion scenarios as it relates to ${\mathrm{\Phi }}_r$. Hence, with the results to be presented, the flow dynamics in L2 of a two-layer stratified medium of known can be described in terms of peak flux and dispersion refraction. As an advantage, our findings could find application in predicting the peak solute plume dispersion when the CBE interface of a pollution control system breaks down (Huo et al., 2022). In addition, incorporating the model into the coding algorithm of reservoir simulators could increase their efficiency and effectiveness in predicting displacement patterns, particularly in smart injection wells management systems for EOR. This could cater for stratification zones of limited hydraulic properties’ knowledge, i.e. zones where only porosity is known (Brouwer et al., 2001; Dawe & Marcelle-de Silva, 2005; Obibuike et al., 2022; Shahrokhi et al., 2014). Other broad application areas include environments where stratifications are often encountered, such as geothermal energy extraction, geological storage of hydrogen and carbon dioxide, chemical separation of mixtures and drug delivery in living organisms (Dejam & Hassanzadeh, 2022).

The subsequent sections of this paper are organised as follows. The Methodology section discusses the data source, the ISE-LHPM setup, the data obtained and the analysis methods applied. After which, the Results and Discussion are presented, followed by the Summary and Conclusions.

2. Methodology

Figure 3 shows the flowchart of the methodology applied in this study. It starts with obtaining the data from the ISEs (section 2) and terminates with summary and conclusions (section 4).

Figure 3. Methodology flowchart.

The data used in this work were obtained from Popoola, Adegoke, Alabi et al. (2010a), Popoola et al. (2010b) and Aubertin et al. (2009). The setup is illustrated in Figure 4. Each of the five LHPMs comprises two layers of different soil samples of known ϕ yielding ${\mathrm{\Phi }}_r$ of 0.8325, 0.8667, 0.9100, 0.9148, and 0.9523, respectively. The experiment setup comprised a 108.5 cm long inlet pipe of radius 2.23 cm and five equal outlet pipes of radii 0.3 cm each. The outlet pipes were joined to the inlet pipe at outlet angles (θ) of 0°, 20°, 50°, 70°, and 90° from the normal. As a control experiment, water was allowed to flow through the empty pipes for 60 s. The discharged water from each outlet was collected separately and used to estimate Q and q, based on Eqn. 5 and 6(6) $\mathit{q}=\frac{Q}{A}=\mathit{K}\cdot \mathit{i}$(6) , respectively (Alabi et al., 2019; H. Wang et al., 2020).

(5) $\mathit{Q}=\frac{V}{60\mathit{s}}$(5)

Figure 4. Schematic diagram of the experiment setup (Popoola et al., 2010a).

where Q is the volumetric flow rate, v is the volume, and 60 s is the amount of time allowed.

(6) $\mathit{q}=\frac{Q}{A}=\mathit{K}\cdot \mathit{i}$(6)

where q is the volume flux, A is the pipe’s cross-sectional area, K is the hydraulic conductivity, i is the hydraulic gradient, and all other notations retain their meaning.

This process was repeated at different inclination angles (α) in order to observe the effect on the flow process. The process described above was done with each of the LHPMs.

2.1. Data

Table 1–4 gives the obtained q at different $\mathit{\alpha }$ ranging from 0 ${ }^{\circ }$ to 20 ${ }^{\circ }$ for ${\mathrm{\Phi }}_r$ 0.8325–0.9524. From the tables, it was observed that the $q_{\mathit{max}}$ are obtained at lower θ (refraction angles from the normal) as the ${\mathrm{\Phi }}_r$ increases.

Table 1. Volume flux at different inclination angles and outlet angles for ${\mathrm{\Phi }}_A=0.8325$

$\mathit{\alpha }$(°)	θ
	0° $\mathit{q}\times {10}^{-3}$ ms^−1	20° $\mathit{q}\times {10}^{-3}$ ms^−1	50° $\mathit{q}\times {10}^{-3}$ ms^−1	70° $\mathit{q}\times {10}^{-3}$ ms^−1	90° $\mathit{q}\times {10}^{-3}$ ms^−1
0	0	0	0	0	0
5	0	0	4.01	1.61	0.43
10	0	0	4.95	3.73	1.85
15	0.90	1.02	5.07	4.91	3.65
20	1.26	2.59	11.43	10.96	4.56

Table 2. Volume flux at different inclination angles and outlet angles for ${\mathrm{\Phi }}_B=0.8667$

α(°)	θ
	0° $\mathit{q}\times {10}^{-3}$ ms^−1	20° $\mathit{q}\times {10}^{-3}$ ms^−1	50° $\mathit{q}\times {10}^{-3}$ ms^−1	70° $\mathit{q}\times {10}^{-3}$ ms^−1	90° $\mathit{q}\times {10}^{-3}$ ms^−1
0	0	0	0	0	0
5	0	6.13	1.89	0.90	0.61
10	0.16	9.10	4.56	2.71	2.04
15	1.34	12.45	7.35	6.09	5.35
20	3.07	15.52	12.65	7.33	6.09

Table 3. Volume flux at different inclination angles and outlet angles for ${\mathrm{\Phi }}_C=0.9100$

$\mathit{\alpha }$(°)	θ
	0° $\mathit{q}\times {10}^{-3}$ ms^−1	20° $\mathit{q}\times {10}^{-3}$ ms^−1	50° $\mathit{q}\times {10}^{-3}$ ms^−1	70° $\mathit{q}\times {10}^{-3}$ ms^−1	90° $\mathit{q}\times {10}^{-3}$ ms^−1
0	0	0	0	0	0
5	0.2	2.71	1.85	1.43	0
10	2.43	5.07	3.14	1.46	0.12
15	6.11	8.4	4.29	2.64	1.02
20	7.19	8.67	4.36	3.92	1.61

Table 4. Volume flux at different inclination angles and outlet angles for ${\mathrm{\Phi }}_D=0.9148$

$\mathit{\alpha }$(°)	θ
	0° $\mathit{q}\times {10}^{-3}$ ms^−1	20° $\mathit{q}\times {10}^{-3}$ ms^−1	50° $\mathit{q}\times {10}^{-3}$ ms^−1	70° ms^−1	90° $\mathit{q}\times {10}^{-3}$ ms^−1
0	0	0	0	0	0
5	2.48	0.9	0.9	0.75	0
10	5.13	1.89	2.43	1.89	0
15	8.31	2.71	2.64	2.43	0.12
20	8.41	7.35	5.11	4.29	0.12

2.2. Data analysis

Table 6 shows the summary data from Tables 1–5 of the ${\theta }_{\mathit{max}}$ at each α, for each of the five LHPMs. Since it was observed that at $\mathit{\alpha }=0^{\circ }$, $\mathit{q}=0$, the whole row was omitted in Table 6.

Table 5. Volume flux at different inclination angles and outlet angles for ${\mathrm{\Phi }}_E=0.9524$

$\mathit{\alpha }$(°)	θ
	0° $\mathit{q}\times {10}^{-3}$ ms^−1	20° $\mathit{q}\times {10}^{-3}$ ms^−1	50° $\mathit{q}\times {10}^{-3}$ ms^−1	70° $\mathit{q}\times {10}^{-3}$ ms^−1	90° $\mathit{q}\times {10}^{-3}$ ms^−1
0	0	0	0	0	0
5	7.88	2.43	1.26	0	0
10	11.02	5.35	4.29	0.2	0
15	16.82	8.41	7.19	2.48	0.12
20	17.8	8.57	8.31	6.01	0.16

Table 6. ${\mathrm{\Phi }}_r$ and their associated ${\theta }_{\mathit{max}}$ at each $\mathit{\alpha }$ level

$\mathit{\alpha }$(°)		${\mathrm{\Phi }}_r$
		0.8325	0.8667	0.9100	0.9148	0.9523
5	${\theta }_{\mathit{max}}$	50°	20°	23°	0°	0°
10		53°	23°	20°	0°	0°
15		60°	23°	18°	0°	0°
20		60°	30°	15°	0°	0°

From Table 6, four plots (Fig. 5–8) and associated models were generated at each α. Furthermore, a combined plot (Fig. A1—shown in appendix A) and model (Eqn. 1A)(1A) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(1A) valid at all the α was generated. The residuals of Eqn. 1A(1A) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(1A) were checked for the satisfaction of Ordinary Least Squares (OLS) regression assumptions using the internally standardised residual (Gray & Woodall, 1994). The internally standardised residual checks for homoskedasticity, as shown in Figure A2.

However, the variance observed in Figure A2 was unsatisfactory, indicating one outlier’s presence. Hence, weighted regression was carried out based on Eqn. 7(7) ${\hat{e}}_{\mathit{i}}={\sum }_i^n{w}_i\left(y_i-{\hat{y}}_{\mathit{i}}\right)$(7) (Romano & Wolf, 2017; Suárez et al., 2017). This was done to prevent the removal of the outlying data point, as such points are known to have undue influence on the regression line, leading to an inaccurate model (Dhakal, 2017).

(7) ${\hat{e}}_{\mathit{i}}={\sum }_i^n{w}_i\left(y_i-{\hat{y}}_{\mathit{i}}\right)$(7)

$w_i=\frac{1}{{\sigma }_{{\hat{e}}_{\mathit{i}}}^2}$ (Davidian & Haaland, 1989)

where ${\hat{e}}_{\mathit{i}}$ is the $i^{\mathit{th}}$ residual, $w_i$ is the $i^{\mathit{th}}$ weight, $y_i$ is the $i^{\mathit{th}}$ the observed value, and ${\hat{y}}_{\mathit{i}}$ is the $i^{\mathit{th}}$ predicted value.

Since only one iteration of the weighted regression yielded improvement in the residual plots, no further iteration was attempted. This resulted in the weighted combined model of Eqn. 12(12) ${\theta }_{\mathit{max}}=3599.66-7579.23{\mathrm{\Phi }}_r+3989.59{{\mathrm{\Phi }}_r}^2$(12) . Figure 9 shows the weighted model superimposed on the scatter plot, while Fig. 10–11 show the internally standardised residual plot and normality plot, respectively.

Fig. 10–11 are only visual tests. Therefore, more formal tests were conducted to ensure that the OLS assumptions were truly satisfied. The tests conducted are the Shapiro–Wilk test for normality (Table 7), at a significance level of 0.05, and the Bartlett test for variance homogeneity at a confidence level of 95% (Devore, 2011; Freund et al., 2010; King & Eckersley, 2019; Odoi et al., 2022). As a precaution, it was ensured that the Shapiro–Wilk test was carried out to verify that the residuals were genuinely normally distributed before the Bartlett test was conducted. The precaution was taken because the Bartlett test only performs optimally on normally distributed data (King & Eckersley, 2019; Odoi et al., 2022).

Table 7. Shapiro–Wilk test (confidence level = 0.05) result summary for the residual

Null Hypothesis:	The sample was drawn from a population that follows a normal distribution
Alternative Hypothesis:	The sample was drawn from a population that does not follow a normal distribution
Sample Size	Computed Statistic	Computed p-value
15	0.8988	0.0920
Result:	Accepted: This statistical test does not provide enough evidence to conclude that the null hypothesis is false.

3. Result and discussion

From Table 6, it can be observed that there is a decreasing trend in ${\theta }_{\mathit{max}}$ with increasing ${\mathrm{\Phi }}_r$, except at $\mathit{\alpha }=5^{\circ }$, with a minor difference in the trend at ${\mathrm{\Phi }}_r=0.9100$. Of course, this is expected since the tendency of the flow path to be disturbed increases with a decrease in porosity viz-a-viz permeability. It can be inferred that with an increasing ${\mathrm{\Phi }}_r$ (${\mathrm{\Phi }}_1$ held constant while ${\mathrm{\Phi }}_2$ reduces, or ${\mathrm{\Phi }}_2$ held constant while ${\mathrm{\Phi }}_1$ increases) ${\theta }_{\mathit{max}}$ will tend to refract towards the normal, i.e., the $q_{\mathit{max}}$ flowline will tend towards the normal. This agrees with the numerical analysis findings of Sebben and Werner (2016a, 2016b) based on K, where they also reported that a less concentrated peak solute plume is more displaced than a more considerably concentrated type. Therefore, in the waterflooding scenario, the flowlines will maintain $q_{\mathit{max}}$ in zones located closer to the contact interface’s normal. At the same time, the peak solute plume flux will undergo a lower dispersal or dilution rate due to the refraction towards the normal (Sebben & Werner, 2016a). In comparison, with a decreasing ${\mathrm{\Phi }}_r$ (held constant while ${\mathrm{\Phi }}_2$ increases, or ${\mathrm{\Phi }}_2$ held constant while ${\mathrm{\Phi }}_1$ decreases) ${\theta }_{\mathit{max}}$ will refract from the normal, i.e., the peak fluid flux flow will refract from the normal. Hence, in a waterflooding scenario, the flowlines will maintain $q_{\mathit{max}}$ in zones located away from the contact interface’s normal as long as the ${\mathrm{\Phi }}_r$ of the layers in contact is low enough. On a similar note, for fluids containing solute plumes (single-phase), e.g., EOR fluids containing nanoparticles or other chemicals, the peak solute plume flux $\left({\theta }_{\mathit{max}}\right)$ will get dispersed towards a higher refraction angle from the normal. This would lead to a faster dilution of the solutes in the L2. If the ${\mathrm{\Phi }}_r$ is low enough, up to 100 times concentration lesser than the concentration in L1 could be observed in L2 (Sebben & Werner, 2016a). These propositions remain valid as long as the flow occurs in an AF system or arrangement.

Furthermore, from Tables 1–6, it is worth noting that there is no definite trend between $\mathit{\alpha }$ and ${\theta }_{\mathit{max}}$. But, $q_{\mathit{max}}$ increases as $\mathit{\alpha }$ increases. These trends indicate that the refraction angle of the peak solute plume, which causes the spread in L2, does not depend on the layer inclination. Therefore, the refraction angle is purely an interfacial property, as reported in the generalisation of Taylor dispersion theory by Dejam and Hassanzadeh (2022). However, the layer inclination does affect the flux rate (volume of the solute plume per unit area per second) through L2. The latter trend is believed to be due to the known sine relationship of $\mathit{a}=\mathit{g}sin\mathit{\alpha }$, where a is the acceleration of the fluid particles, g is the acceleration due to gravity, and $\mathit{\alpha }$ retains its meaning (Ceccarelli et al., 2018).

3.1. Individual models

From Figure 5, at $\mathit{\alpha }=5^{\circ }$, it is observed that with increasing ${\mathrm{\Phi }}_r$, ${\theta }_{\mathit{max}}$ decreases. The fitted model is:

(8) ${\theta }_{\mathit{max}}=2372.19-4897.88{\mathrm{\Phi }}_r+2528.95{{\mathrm{\Phi }}_r}^2$(8)

where ${\theta }_{\mathit{max}}$ is the refraction angle of the peak solute plume concentration, and ${\mathrm{\Phi }}_r$ is the porosity ratio of the layers in contact.

Eqn. 8(8) ${\theta }_{\mathit{max}}=2372.19-4897.88{\mathrm{\Phi }}_r+2528.95{{\mathrm{\Phi }}_r}^2$(8) has a coefficient of determination (R²) of 0.82. It implies that 82% of the dependent variable $\left({\theta }_{\mathit{max}}\right)$ can be explained by the independent variable $\left({\mathrm{\Phi }}_r\right)$.

Figure 5. ${\theta }_{\mathit{max}}$ against ${\mathrm{\Phi }}_r$ at $\mathit{\alpha }=5^{\circ }$.

From Figure 6, at $\mathit{\alpha }={10}^{\circ }$, it is observed that with increasing ${\mathrm{\Phi }}_r$, ${\theta }_{\mathit{max}}$ decreases. The fitted model is:

(9) ${\theta }_{\mathit{max}}=2798.60-5818.35{\mathrm{\Phi }}_r+3024.78{{\mathrm{\Phi }}_r}^2$(9)

where all notations retain their meaning.

Eqn. 9ʹs $R^2=0.89$. It can then be said that 89% of ${\theta }_{\mathit{max}}$ can be explained by ${\mathrm{\Phi }}_r$, at $\mathit{\alpha }={10}^{\circ }$

Figure 6. ${\theta }_{\mathit{max}}$ against ${\mathrm{\Phi }}_r$ $\mathit{\alpha }={10}^{\circ }$.

As observed in Figures 5–6, in Fig. 7, ${\theta }_{\mathit{max}}$ decreases with increasing ${\mathrm{\Phi }}_r$. However, a trend reversal is observed at about ${\mathrm{\Phi }}_r=0.94$. Also, compared to the decrease in ${\theta }_{\mathit{max}}$ observed in Figures 5–6, the decrease observed in Fig. 7 is steeper. The model generated is:

(10) ${\theta }_{\mathit{max}}=3845.96-8117.14{\mathrm{\Phi }}_r+4284.29{{\mathrm{\Phi }}_r}^2$(10)

where all notations retain their meaning.

Eqn. 10ʹs $R^2=0.92$. It is of better fit based on the $R^2$ value when compared to those obtained in Figures 5–6, since 92% of the ${\theta }_{\mathit{max}}$ can be explained by the ${\mathrm{\Phi }}_r$ data.

Figure 7. ${\theta }_{\mathit{max}}$ against ${\mathrm{\Phi }}_r$ at $\mathit{\alpha }={15}^{\circ }$

From Figure 8, Eqn. 11(11) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(11) was obtained as the model:

(11) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(11)

where all notations retain their meaning.

Eqn. 11ʹs $R^2=0.96$. It is the best model based on $R^2$ value

Figure 8. ${\theta }_{\mathit{max}}$ against ${\mathrm{\Phi }}_r$ at $\mathit{\alpha }={20}^{\circ }$.

As observed from Eqns. 8(8) ${\theta }_{\mathit{max}}=2372.19-4897.88{\mathrm{\Phi }}_r+2528.95{{\mathrm{\Phi }}_r}^2$(8) –11(11) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(11) , the lowest R² being approximately 0.82, all the models are of good fit. However, R² can be biased in small sample size conditions like we have (Akossou & Palm, 2013). Nevertheless, in comparison to other statistical alternatives like the Index of Agreement (d) and the Maximal Information Coefficient (MIC), it provides a better result (Valbuena et al., 2019). In addition, to have a reasonable reliance on the models, it is necessary to verify the suitability of the models’ form for our context of use (Non-Darcian flow; Howarth, 2017; Oloro, 2020; Valbuena et al., 2019). Deng and Li (2021) and Hansbo (2001, 1960) agree with their suitability and conclude that second-order polynomials have smaller Mean Square Errors (MSE) and a high R² with experimental data.

3.2. Combined model

From Figure 9, a gradual decrease in ${\theta }_{\mathit{max}}$ with increase in ${\mathrm{\Phi }}_r$ is observed. It follows the same trend as the individual plots (Figures 5–8). The superimposed weighted model line is:

(12) ${\theta }_{\mathit{max}}=3599.66-7579.23{\mathrm{\Phi }}_r+3989.59{{\mathrm{\Phi }}_r}^2$(12)

where all notations retain their meaning.

Eqn. 12ʹs $R^2=0.99$. The $R^2$ value $\left(R^2=0.99\right)$ for the weighted model (Eqn. 12)(12) ${\theta }_{\mathit{max}}=3599.66-7579.23{\mathrm{\Phi }}_r+3989.59{{\mathrm{\Phi }}_r}^2$(12) is an improvement over the unweighted model’s (Eqn. 1A(1A) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(1A) ) $R^2$ value $\left(R^2=0.89\right)$. Also, Eqn. 12ʹs internally standardised residuals plot (Figure 10) shows improved residual variance homoscedasticity compared to Eqn. 1A’s (Figure A2). As seen in Figure 10, there are no outliers (data points beyond −2 or 2), which indicates the absence of any data point having undue influence on the regression line (Peck et al., 2008)

Figure 9. ${\theta }_{\mathit{max}}$ against ${\mathrm{\Phi }}_r$ at $\mathit{\alpha }=0^{\circ }-{20}^{\circ }$.

Figure 10. Internally standardised residual plot for the weighted combined model.

Figure 11. Residual normal probability plot for the weighted combined model.

Figure 11 shows the normal probability plot of the residuals. It indicated some divergence from normality. However, from Table 7, the Shapiro–Wilk test statistically confirmed the residuals as normally distributed. From the Bartlett test for homogeneity, a statistic of 2.95 was obtained, giving an associated P-value of 0.229. The Bartlett test result indicates that there is insufficient statistical evidence to reject the null hypothesis of the residuals being homoscedastic at the chosen P-value of 0.05. Hence, based on the satisfaction of OLS assumptions, the weighted model (Eqn. 12)(12) ${\theta }_{\mathit{max}}=3599.66-7579.23{\mathrm{\Phi }}_r+3989.59{{\mathrm{\Phi }}_r}^2$(12) is efficient and the Best Linear Unbiased Estimator (BLUE; Pötscher & Preinerstorfer, 2022).

4. Summary and conclusions

A study has been carried out to highlight the flow dynamics as it relates to the refraction angle $\left({\theta }_{\mathit{max}}\right)$ or dispersion of the peak volume or solute plume flux of a single-phase fluid flowing through the contact interface of two layers of Layered Heterogeneous Porous Media (LHPM) characterised by different porosities. The LHPM is an idealised representation of a stratified medium. The flow was orientated so that it originates from a lower porosity vis-à-vis permeability layer (top layer—L1) and flows into a higher porosity vis-à-vis permeability layer (bottom layer—L2).

From the analysis of the data, it can be inferred that:

1. Flowlines will maintain $q_{\mathit{max}}$ (peak volume or solute plume flux) in zones located away from the contact interface’s normal, as long as the porosity ratio $\left({\mathrm{\Phi }}_r\right)$ of the layers that make up such a contact interface is low enough. This would lead to a higher dispersal angle and fast rate of dilution of solute plumes in the bottom layer. On the other hand, with ${\mathrm{\Phi }}_r$ tending to unity, $q_{\mathit{max}}$ flowline will tend towards the normal. In this case, a lower dispersal angle and slow rate of dilution of solute plumes in the bottom layer would be observed.

2. The inclination angle $\left(\mathit{\alpha }\right)$ of the stratification units or layers does not correlate with ${\theta }_{\mathit{max}}$. However, $\mathit{\alpha }$ correlates with $q_{\mathit{max}}$

3. For an AF system, the best-unbiased estimator model that could be used to predict ${\theta }_{\mathit{max}}$ based on ${\mathrm{\Phi }}_r$ is:

${\theta }_{\mathit{max}}=3599.66-7579.23{\mathrm{\Phi }}_r+3989.59{{\mathrm{\Phi }}_r}^2$

where ${\theta }_{\mathit{max}}$ is the refraction angle of the $q_{\mathit{max}}$ or the peak solute plume concentration and ${\mathrm{\Phi }}_r$ is the porosity ratio of the contact interface.

• (4)The results of this work could find application in any porous environment characterised by stratification with an AF system flow regime. Such instances include the mapping of ${\theta }_{\mathit{max}}$ and probable dilution rate of flooded fluids in optimised smart injection well management systems for EOR. In addition, the model could optimise CBE covers used for groundwater protection by predicting the path of the peak solute flux, which could then be considered to improve the covers’ performance.

This study was based on a single-phase fluid undergoing a one-way directional flow in an AF system. However, in most cases, such dynamics do not restrain flow in nature. Occurrence of the reverse flow, i.e. DF, is possible. Flow could also be in multi-phase like those experienced in hydrocarbon reservoirs. Hence, these limitations present possible future research opportunities.

Nomenclature

${\mathrm{\Phi }}_r$ Porosity ratio

K Hydraulic conductivity

V Volume of fluid

i Hydraulic gradient

$q_{\mathit{max}}$ Peak or maximum volume flux of solute plume

q Volume flux (m/s)

Q Volumetric flow rate (m³s⁻¹)

A Cross-sectional area perpendicular to the line offlow (m²)

L1 Layer 1 or the top layer for the AF system whereflow originates

L2 Layer 2 or the bottom layer for the AF system where flow terminates

*L1 Layer 1 or the top layer for the DF system where flow originate

*L2 Layer 2 or the bottom layer for the DF system where flow terminates

Greek symbols

α Inclination angle

θ Flowline deflection angle from the normal (Outlet angle) or refraction angle of solute plume flux

${\theta }_1$ Entry angle

${\theta }_2$ Refraction angle

${\theta }_{\mathit{max}}$ Refraction angle from the normal of the flowline associated with or refraction angle of peak solute plume flux

${\mathrm{\Phi }}_1$ L1 porosity

${\mathrm{\Phi }}_2$ L2 porosity

$\ast {\mathrm{\Phi }}_1$*L1 porosity

$\ast {\mathrm{\Phi }}_2$ *L2 porosity

ϕ Porosity

${\mathrm{\Phi }}_r$ Porosity ratio

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The authors received no direct funding for this research.

Appendix A

The superimposed model line is:

(1A) ${\theta }_{\mathit{max}}=3601.24-7536.03{\mathrm{\Phi }}_r+3942.25{{\mathrm{\Phi }}_r}^2$(1A)

where all notations retain their meaning.

Eqn. 1A’s $R^2=0.89$

Figure A1. ${\theta }_{\mathit{max}}$ against ${\mathrm{\Phi }}_r$ at $\mathit{\alpha }=0^{\circ }-{20}^{\circ }$

Figure A2. Internally standardised residual plot for the unweighted combined model.