4.1. Parameter selection

The selection of parameter values in this study is based on a comprehensive review of existing literature (Chen & Wu Reference Chen and Wu2012). Table 1 presents the range of parameter values used in these investigations. The Darcy number ( $ D_a$ ) indicates permeability, where a higher $ D_a$ signifies more permeable, less compact pore geometry, and a lower $ D_a$ signifies less permeable, more compact geometry. The viscosity ratio ( $ M$ ), representing the ratio of two viscosities, is crucial for exploring the impact on fluid dynamics. The pressure gradient ( $ Q$ ), expressed as $ -({{\rm d}p}/{{\rm d}x)}/({\mu U}/{H^2})$ , allows for the analysis of the influence on the composite region, with different $ Q$ values representing varying pressure gradients, facilitating the examination of flow behaviour under different pressure conditions. Finally, time ( $ t^*$ ) represents the temporal aspect of pollutant dispersion, enabling the analysis of how pollutant concentration and distribution evolve over time. The Péclet numbers ( ${\textit{Pe}}_p$ and ${\textit{Pe}}_{\!f}$ ) are dimensionless numbers that characterise the relative importance of advection to diffusion in the transport processes, with their values indicating the dominance of one mechanism over the other.

Table 1.

Range of parameter values.

4.2. Computational results and validations

The expression for the two-dimensional concentration distribution can be derived from (3.17) and (3.18) in the new coordinate system $\{\zeta ,y\}$ as

(4.1) \begin{equation} C=\begin{cases} C^{(0)}+\epsilon {\textit{Pe}}_{p}N_p(y)\dfrac {\partial C^{(0)}}{\partial \zeta }+\epsilon ^2{\textit{Pe}}_p X_p(y)\dfrac {\partial ^2 C^{(0)}}{\partial {\zeta }^2}, & 0\lt y\leqslant h,\\[10pt] C^{(0)}+\epsilon {\textit{Pe}}_{\!f}N_{\!f}(y)\dfrac {\partial C^{(0)}}{\partial \zeta }+\epsilon ^2{\textit{Pe}}_{\!f} X_{\!f}(y)\dfrac {\partial ^2 C^{(0)}}{\partial {\zeta }^2}, & h \leqslant y\lt 1. \end{cases} \end{equation}

We used the finite difference discretisation technique, as explained in Appendix D, to obtain the following discrete forms of (3.7) and (3.8) as follows:

(4.2) \begin{align} & {\textit{Pe}}_p \left \langle \frac {1}{\textit{Pe}} \right \rangle \frac {C_{p_{i,j}}^{n+1} - C_{p_{i,j}}^{n}}{\Delta t^*} \nonumber\\ & \quad + {\textit{Pe}}_p \left \{ {\text{max}} (u_p - \langle u \rangle , 0) \frac {C_{p_{i,j}}^{n} - C_{p_{i-1,j}}^{n}}{\Delta \zeta } - {\text{max}} (\langle u \rangle - u_p, 0) \frac {C_{p_{i+1,j}}^{n} - C_{p_{i,j}}^{n}}{\Delta \zeta } \right \} \nonumber\\ & = \frac {C_{p_{i+1,j}}^{n} - 2C_{p_{i,j}}^{n} + C_{p_{i-1,j}}^{n}}{\Delta \zeta ^2}+\frac {C_{p_{i,j+1}}^{n} - 2C_{p_{i,j}}^{n} + C_{p_{i,j-1}}^{n}}{\Delta y^2}, \quad j=0,1,\ldots ,j_h, \end{align}

(4.3) \begin{align} & {\textit{Pe}}_{\!f} \left \langle \frac {1}{\textit{Pe}} \right \rangle \frac {C_{f_{i,j}}^{n+1} - C_{f_{i,j}}^{n}}{\Delta t^*} \nonumber\\ & \quad + {\textit{Pe}}_{\!f} \left \{ {\text{max}} (u_p - \langle u \rangle , 0) \frac {C_{f_{i,j}}^{n} - C_{f_{i-1,j}}^{n}}{\Delta \zeta } - {\text{max}} (\langle u \rangle - u_{\!f}, 0) \frac {C_{f_{i+1,j}}^{n} - C_{f_{i,j}}^{n}}{\Delta \zeta } \right \} \nonumber\\ & = \frac {C_{f_{i+1,j}}^{n} - 2C_{f_{i,j}}^{n} + C_{f_{i-1,j}}^{n}}{\Delta \zeta ^2}+\frac {C_{f_{i,j+1}}^{n} - 2C_{f_{i,j}}^{n} + C_{f_{i,j-1}}^{n}}{\Delta y^2}, \quad j=j_h,j_h+1,\ldots N_y. \end{align}

Here $i=0,1,2,\ldots ,M_x$ and $j=0,1,2\ldots ,N_y$ ; $M_x$ and $N_y$ represent the maximum number of grid points along the $\zeta$ and $y$ direction, respectively. The computational domain considered is $[-20,20]\times [0,1]$ . We discretise the domain using $2001$ grid points in the longitudinal direction and $101$ grid points in the transverse direction, corresponding to $\Delta x = L_x/(M_x-1)$ and $\Delta y = L_y/(N_y-1)$ . The time step is chosen as

(4.4) \begin{equation} \Delta t = \frac {0.4}{(1/\Delta x^2 + 1/\Delta y^2)}, \end{equation}

which satisfies the stability condition of the explicit scheme.

Figure 2.

Comparison between analytical and numerical results of a two-dimensional concentration distribution for different times with $M=1$ , ${\textit{Pe}}_p={\textit{Pe}}_{\!f}=1$ , $Q=1$ , $h=0.5$ , $D_a=0.1$ .

In figure 2 we present a comparative analysis between the analytical and numerical results over time. The graph shows a strong relation between the two-dimensional concentration distributions obtained from both the analytical and computational approaches. Although experimental data for this specific case is not available, we have provided a comparison with results derived from computational methods. This comparison highlights the consistency between the two methods and validates the reliability of the computational model in capturing the behaviour predicted by the analytical solution. Despite the absence of experimental verification, the close agreement between these solution methods reinforces the accuracy and applicability of both methods in predicting this phenomenon.

4.3. Dispersion coefficients

The effective dispersion coefficient can be expressed in (3.57) as

(4.5) \begin{equation} D_T = 1 + D_T^*. \end{equation}

The value of $D_T^*$ takes into account the impact of lateral diffusion and advection and is calculated as $-{\langle PeuN(y) \rangle }/{\langle {1}/{\textit{Pe}} \rangle }$ . On the right-hand side, the constant term of (4.5), i.e. the one representing the effect of longitudinal diffusion, can be safely ignored to focus on the enhanced effects of lateral diffusion and advection. The dispersion coefficient ( $D_T^*$ ) is influenced by various parameters such as ${\textit{Pe}}_{\!f}$ , ${\textit{Pe}}_p$ , $M$ , $h$ , $Q$ and $D_a$ . When the Darcy number( $D_a$ ) is large, $D_T^*$ becomes 0.0083 $Pe^2_{\!f}$ , which is the same as the finding of Barik & Dalal (Reference Barik and Dalal2018). This value represents the degree of dispersion of the solute in the fluid medium under single-zone open channel conditions (i.e. $h=0$ in our case).

Figure 3.

Variation of dispersion for (a) $M=1, {\textit{Pe}}_p={\textit{Pe}}_{\!f}=1, Q=1$ and (b) $D_a=1, {\textit{Pe}}_p=1, M=1, Q=1$ .

Figure 3(a) shows the variation of dispersion for different $D_a$ values. We can see from this that an increase in $D_a$ leads to decrease in the dispersion value and a flatter dispersion graph, which is same as in Chen & Wu (Reference Chen and Wu2012). When the value of $D_a$ is 0.01, the dispersion coefficient displays a discernible pattern. It gradually increases until it reaches a specific value of $h$ , which is not fixed for every $D_a$ , after which it experiences a rapid decrease. When $D_a$ is high, fluid flow through porous media becomes more efficient, and advection becomes increasingly essential compared with diffusion. When advection dominates (high $D_a$ ), the dispersion graph becomes flatter because advection rapidly homogenises the solute concentration within the porous medium. The variation of dispersion graphs with different Péclet numbers for the fluid phase ( ${\textit{Pe}}_{\!f}$ ), as depicted in figure 3(b), illustrates a clear relationship between ${\textit{Pe}}_{\!f}$ values and the dispersion coefficient. As ${\textit{Pe}}_{\!f}$ increases, it can be noticed that the dispersion coefficient also increases (Chen & Wu Reference Chen and Wu2012).

Figure 4(a) illustrates the variation of dispersion coefficient with respect to the height of the porous region, $h$ , for different values of $M$ , which represents the ratio of effective viscosity to the dynamic viscosity. For each value of $M$ , the dispersion initially increases as $h$ increases from 0 to 0.3, and then shows varying behaviour for higher $h$ . When $h=0.5$ , the dispersion remains relatively stable or slightly decreases, but for higher values of $h$ (around 0.8 and 1), we observe a general decline in dispersion. The height of the porous region, $h$ , affects the flow dynamics and the spread of particles within the porous medium. At lower heights ( $h = 0$ to $h = 0.3$ ), the flow is less constrained because the porous region is not fully developed, leading to weaker interactions between the fluid and porous medium. As a result, the dispersion initially increases due to less resistance and more mobility of fluid particles. However, as $h$ increases, the porous medium becomes more dominant, imposing greater resistance to flow, which in turn slows down the dispersion. This leads to a peak or plateau in dispersion values around $h = 0.3$ , followed by a decline as the resistance further increases with larger $h$ . This phenomenon can be seen in the values displayed in table 2.

Table 2.

Values of $D_T^*$ for different values of $M$ and $h$ .

Figure 4.

Variation of dispersion for (a) $ {\textit{Pe}}_p={\textit{Pe}}_{\!f}=1, Q=1$ and $D_a=1$ , and (b) $ {\textit{Pe}}_p={\textit{Pe}}_{\!f}=1, M=1$ and $D_a=1$ .

Table 3 and figure 4(b) present the variations in dispersion coefficient ( $D_T^*$ ) as influenced by the height ( $ h$ ) of the porous region and the pressure gradient ( $ Q$ ). With a constant pressure gradient, the dispersion coefficient increases with the increase in the height of the porous region, and with a fixed height of the porous medium, the dispersion coefficient increases as the pressure gradient increases.

Table 3.

Dispersion ( $D_T^*$ ) values at various heights ( $h$ ) for different pressure gradients ( $Q$ ).

4.4. Mean concentration

The mean concentration distribution can be computed using (3.17), (3.18) and (3.72) as

(4.6) \begin{equation} \langle C\rangle =C^{(0)} = \frac {1}{\sqrt { 4\pi D_T t^*}}\exp \left (\frac {-\zeta ^2}{4 D_T t^*}\right )\!. \end{equation}

Figure 5.

Longitudinal distribution of mean concentration distributions for different values of (a) time, where $M=1$ , ${\textit{Pe}}_p={\textit{Pe}}_{\!f}=1$ , $Q=1$ , $D_a=0.1$ and (b) Péclet numbers ratio, where $t^*=1$ , $Q=1$ , $M=1$ , $D_a=0.1$ .

Figure 6.

Longitudinal distribution of mean concentration distributions for $M=1,{\textit{Pe}}_P={\textit{Pe}}_{\!f}=1,Q=1$ and $t^*=1$ .

The variations of the longitudinal distribution of mean concentration can be seen from figures 5, 6 and 7. Figure 5(a) displays the change in the mean concentration distribution with time. At a fixed time, the mean concentration exhibits a Gaussian distribution, which is characterised by its symmetrical shape and adherence to statistical normality. The axis of symmetry remains unchanged, therefore partitioning the distribution into two symmetrical components. Over time, the apex of the distribution progressively decreases and becomes flatter.

The graph in figure 5(b) demonstrates how the ratio of the Péclet numbers $({\textit{Pe}}_{\!f}/{\textit{Pe}}_p)$ affects the mean concentration at a specific height $(h = 0.5)$ . When $({\textit{Pe}}_{\!f}/{\textit{Pe}}_p)$ increases, the Péclet number of the fluid phase $({\textit{Pe}}_{\!f})$ becomes higher than that of the porous phase $({\textit{Pe}}_p)$ , resulting in a decrease in the peak of the mean concentration graph. This occurs because advection in the fluid phase becomes more influential, resulting into better transport of solutes. As a result, the solutes are more effectively carried away from the peak concentration point, causing a lower peak.

Figure 7.

Longitudinal distribution of mean concentration distributions for $M=1,{\textit{Pe}}_p={\textit{Pe}}_{\!f}=1,h=0.5,D_a=0.1$ and $t^*=1$ .

Figures 6(a) and 6(b) show that the mean concentration distribution undergoes significant changes with the Darcy number $(D_a)$ in a two-layered channel flow. It is worth noting that these trends correlate with the behaviour of the dispersion coefficient, as demonstrated in figure 3. In figure 3(a) for $h = 0.3$ , we can see that dispersion increases with a lower Darcy number, causing the peak of the mean concentration to decrease, as shown in figure 6(a). This happens because as $D_a$ goes down, the effect of molecular diffusion on dispersion becomes stronger. This makes solutes spread out more and lowers their peak concentrations. Conversely, for $h = 0.8$ , a lower $D_a$ corresponds to decreased dispersion, resulting in an increased peak of the mean concentration, as illustrated in figure 6(b). In this case, advection is more prominent in the channel, carrying solutes away more quickly and thus reducing peak concentrations. Figure 7 displays the distribution of mean concentration for different $Q$ values while keeping the parameters constant at $h = 0.5$ . Our observation indicates that lower $Q$ values result in a higher concentration peak, thus leading to greater solute spreading when the $Q$ value is high. This finding is consistent with the dispersion graph presented in figure 4(b).

4.5. Longitudinal concentration distribution

The expression for the two-dimensional concentration distribution up to second order can be derived from (3.17) and (3.18) in the new coordinate system $\{\zeta ,y\}$ as

(4.7) \begin{equation} C=\begin{cases} C^{(0)}+\epsilon {\textit{Pe}}_{p}N_p(y)\dfrac {\partial C^{(0)}}{\partial \zeta }+\epsilon ^2{\textit{Pe}}_p X_p(y)\dfrac {\partial ^2 C^{(0)}}{\partial {\zeta }^2}, & 0\lt y\leqslant h,\\[8pt] C^{(0)}+\epsilon {\textit{Pe}}_{\!f}N_{\!f}(y)\dfrac {\partial C^{(0)}}{\partial \zeta }+\epsilon ^2{\textit{Pe}}_{\!f} X_{\!f}(y)\dfrac {\partial ^2 C^{(0)}}{\partial {\zeta }^2}, & h \leqslant y\lt 1. \end{cases} \end{equation}

Figures 8 and 9 show the longitudinal solute distribution along the channel, illustrating the gradual transition from zeroth order to second order at different times. Initially, for $t^* \lt 1$ , the contributions of each order term are significant in shaping the solute concentration along the channel. However, as time progresses, the influence of higher-order terms diminishes. After $t^* = 1$ , neglecting terms beyond the second order introduces minimal error, as indicated by figures 8 and 9. As a result, the graphs for zeroth-order, first-order and second-order approximations become nearly indistinguishable for a considerable period when $t^* \gt 1$ .

Figure 8.

Longitudinal distribution of concentration calculated as asymptotic expansions from zeroth order to second order at different times, where $M=1$ , ${\textit{Pe}}_p={\textit{Pe}}_{\!f}=1$ , $h=0.3$ , $Q=1$ , $y=0.5$ , $D_a=0.1$ .

Figure 9.

Longitudinal distribution of concentration calculated as asymptotic expansions from zeroth order to second order at different times, where $M=1$ , ${\textit{Pe}}_p={\textit{Pe}}_{\!f}=1$ , $h=0.8$ , $Q=1$ , $y=0.5$ , $D_a=0.1$ .

Figure 10.

Concentration contour: the horizontal coordinate represents the variable $\zeta$ , while the vertical coordinate represents the variable $y$ for different values of time ( $t^*$ ) and ${{\textit{Pe}}_{\!f}}/{{\textit{Pe}}_p}$ where (a) $M=1, Q=1, D_a=0.1, {\textit{Pe}}_p={\textit{Pe}}_{\!f}=1$ and (b) $M=1, D_a=0.1, t^*=1, Q=1$ .

Figures 10 and 11 display the concentration contours for different values of the parameters. Figure 10(a) presents the solute spreading dynamics over time in the two-layered channel system. Here, we plot the contour graphs for the porous region, which spans from 0 to 0.5 and the fluid region from 0.5 to 1. Initially, the solute is injected across the line $\zeta =0$ , the solute spreads due to advection and diffusion. The solute concentration experiences differential transport rates in the two regions: it travels faster in the fluid region compared with the porous area. This difference in transport velocities is primarily due to the higher permeability and reduced resistance to flow in the fluid region, allowing solute molecules to move more freely. As time progresses, the concentration centroid experiences a reduction in its value, indicating the spreading of the solute throughout the channel. This phenomenon arises from advection, which propels solute downstream, and diffusion, which causes dispersion by random molecular motion, ultimately leading to a broader solute concentration distribution.

In figure 10(b) the observation of a ratio between two Péclet numbers ( ${{\textit{Pe}}_p}/{{\textit{Pe}}_{\!f}}$ ) equals 0.5, resulting in the concentration cloud being at the origin near the bottom plate and near the upper plate along the $y$ axis. As the ratio of these Péclet numbers increases, solute dispersion becomes more pronounced. Furthermore, the observation that the solute cloud spreads more in the fluid region than in the porous region suggests that the fluid layer is more efficient in transporting solutes over long distances due to higher advection effects. These trends align with the behaviour of the dispersion coefficient, which is influenced by the relative contributions of advection and diffusion in the two regions, impacting the overall spread of the solute cloud.

Figure 11.

Concentration contours: the horizontal coordinate represents $\zeta$ , while the vertical coordinate represents $y$ for different values of $D_a$ where (a) $M=1, Q=1, {\textit{Pe}}_p={\textit{Pe}}_{\!f}=1, h=0.3$ and (b) $M=1, Q=1, {\textit{Pe}}_p={\textit{Pe}}_{\!f}=1, h=0.8$ .

The observations in figures 11(a) and 11(b) provide valuable information into the behaviour of solute concentration distribution in a two-layered channel flow with differing heights ( $h = 0.3$ and $h = 0.8$ ). In figure 11(a) at $h = 0.3$ , it is clear that with decreasing Darcy number, solute spreads more. In contrast, at $h = 0.8$ (figure 11 b), the opposite trend is observed, an increase in the Darcy number results in more solute spreading. These observations align with the behaviour of the dispersion coefficient, as a lower Darcy number increases dispersion at $h = 0.3$ , while a higher Darcy number does so at $h = 0.8$ . The consistency between the dispersion and mean concentration results underline the interplay between advection and diffusion in influencing solute transport characteristics in two-layered channel flows with different heights.

4.6. Application of classical Darcy model in porous region

So far we were carrying out our investigations based on the extended Darcy equation (for the flow in porous media). However, based on one of the reviewers’ suggestions, we have now carried out the studies based on pure Darcy equations for the flow with low Darcy numbers and the extended Darcy equation for the flows with greater than those Darcy numbers. Fluid flow through porous media has been characterised using the empirical relationship known as Darcy’s law, i.e.

(4.8) \begin{equation} -\dfrac {{\rm d}p}{{\rm d}x'}-\mu \dfrac {u^{\prime}_{\!p}(y')}{\kappa }=0, \quad 0 \leqslant y'\leqslant h', \end{equation}

and the fluid motion in the fluid region (which is the same as (2.2))

(4.9) \begin{equation} - \dfrac {{\rm d}p}{{\rm d}x'} +\mu \dfrac {{\rm d}^2u^{\prime}_{\!f}}{{\rm d}y'^2} = 0,\quad h' \leqslant y' \leqslant H. \end{equation}

The initial and boundary conditions are

(4.10) \begin{align}&\qquad u^{\prime}_{\!p}(y')=u'_b \quad \text{at} \,y'=h', \end{align}

(4.11) \begin{align}& \dfrac {{\rm d}u^{\prime}_{\!f}(y')}{{\rm d}y'}=\beta '(u'_b-u^{\prime}_{\!p}) \quad \text{at}\,y'=h'. \end{align}

Figure 12.

Dispersion graph of porous media using Darcy’s law for (a) different $Q$ and (b) different ${\textit{Pe}}_{\!f}$ .

Figure 13.

Longitudinal distribution of mean concentration distributions for different values of (a) $h=0.3$ , $t^*=1$ , $D_a=0.001$ , $Q=1$ , $M=1$ and (b) $h=0.8$ , $t^*=1$ , $D_a=0.001$ , $Q=1$ , $M=1$ .

Here also $U$ is the flow velocity of the free surface:

(4.12) \begin{equation} u^{\prime}_{\!f}(y')= U \quad \text{at}\, y'=H. \end{equation}

Solving non-dimensional equations, we get the velocity components in the fluid and porous region, respectively, as

(4.13) \begin{align}& u_{\!f}=-Q\frac {y^2}{2}+ (Qh+\beta (u_b-u_p) )y+\left (1+\frac {Q}{2}-Qh-\beta (u_b-u_p)\right )\!, \end{align}

(4.14) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad u_p=\textit{QD}_a, \end{align}

(4.15) \begin{align}&\qquad\qquad\quad u_b=\frac {Q\dfrac {h^2}{2}-Qh+1+\dfrac {Q}{2}-(h-1)\beta {\textit{QD}}_a}{1-(h-1)\beta }, \end{align}

where $u_b$ is the slip velocity and $\beta =({\alpha }/{\sqrt {\kappa }})$ is a parameter that depends on permeable surface properties; $\alpha$ is a dimensionless parameter that contains structural characteristic of the porous media known as the slip coefficient, defined as $\alpha =\sqrt {({\mu }/{\mu _e})}.$

Figure 14.

Longitudinal distribution of mean concentration distributions for different values of (a) $h=0.3$ , $t^*=1$ , $D_a=0.1$ , $Q=1$ , $M=1$ and (b) $h=0.1$ , $t^*=1$ , $D_a=0.001$ , $Q=1$ , $M=1$ .

Figure 15.

Darcy versus extended Darcy, where $h=0.3$ , $M=1.05$ , $Q=1$ , $t^*=1$ .

After obtaining the velocity profile, the multi-scale method is applied to determine the effective dispersion coefficient $ D^*_T$ and the mean concentration distribution $\langle C \rangle$ . Figures 13 and 14 indicate that when the height $( h )$ is small and permeability is low, the classical Darcy model and the extended Darcy model yield coinciding results. However, as $( h )$ increases or permeability becomes significant, noticeable differences between the two models emerge. From this observation, we can conclude that in cases where the porous medium is very thin and exhibits high permeability, the classical Darcy model is appropriate. Conversely, in scenarios where the porous region is relatively thicker and the permeability is not significantly high, the extended Darcy model (Brinkman model) should be employed for more accurate predictions, supporting the results of Winter et al. (Reference Winter, Valsamidou, Class and Flemisch2022). We present a detailed comparison between the classical Darcy law and the extended Darcy (Brinkman) model by varying the channel thickness parameter $h$ (considered at values $h = 0.3, 0.5, 0.8$ ) and the viscosity ratio $M$ (taken as $M = 1$ and $M = 1.05$ ). For the case $M = 1$ with permeability $\kappa = 10^{-5}$ , the difference between the Darcy and Brinkman models is not clearly visible (figure 13 b). To highlight the distinction more effectively, we next consider the case $M = 1.05$ with $\kappa = 0.1$ (figure 15). For $D_a = 1$ , the predictions of the extended Darcy (Brinkman) model differ markedly from the classical Darcy model, even for thin porous layers ( $h = 0.3$ ). This is because the viscous penetration length is $\mathcal{O}(H)$ , allowing interface shear to penetrate through the entire porous region. The classical Darcy model, lacking this shear-diffusion term, overconfines solute near the core. This aligns with studies showing that the Brinkman term becomes increasingly significant at higher permeability and across boundaries. For $D_a = 0.1$ and $0.01$ , the difference between the two models remains clear. The Brinkman solution continues to exhibit a smoother gradient at the interface and a broader profile. Here, the viscous penetration length $\ell _v$ is smaller but still larger than the porous thickness $h$ , so shear penetration affects the full porous layer. The Darcy model, lacking this mechanism, predicts a sharper decay of concentration towards the wall. When $D_a \leqslant 0.001$ , permeability is low and the drag term dominates. In the Darcy model, the porous layer becomes nearly stagnant, and the concentration drops sharply near the wall. The extended Darcy model still shows slightly higher near-wall concentrations due to limited shear penetration, but because $h$ is small the overall differences shrink compared with the high- $D_a$ case.

In figures 16 and 17, when $D_a = 1$ , the permeability $\kappa$ is such that the viscous penetration length $\ell _v \sim \sqrt {\kappa }$ is of order $H$ . This allows shear stresses from the free-fluid region to travel deep into the porous layer. The extended Darcy (Brinkman) model, which includes the $\mu _e{\nabla} ^2\boldsymbol{u}$ viscous diffusion term, transmits momentum into most of the porous thickness, sustaining a small but finite velocity there and enhancing cross-layer advection.The classical Darcy model neglects viscous shear inside the porous medium, producing a rapid velocity drop near the interface and confining solute transport largely to the free-fluid core. Because $h$ is larger than in figure 15, the effect of Brinkman shear penetration is more pronounced, yielding broader $\langle C \rangle (\zeta )$ profiles and a lower centreline peak.

Figure 16.

Darcy versus extended Darcy, where $h=0.5$ , $M=1.05$ , $Q=1$ , $t^*=1$ .

Figure 17.

Darcy versus extended Darcy, where $h=0.8$ , $M=1.05$ , $Q=1$ , $t^*=1$ .

For $D_a = 0.1$ and $0.01$ , the viscous penetration length $\ell _v$ decreases but remains greater than $h$ , allowing shear penetration to influence the entire porous layer. The Brinkman model still produces smoother gradients at the fluid–porous interface and broader concentration distributions than the Darcy model.

At low $D_a$ , permeability is small and the drag term dominates. In the Darcy model the porous layer becomes nearly stagnant, resulting in a steep drop in concentration near the wall. The extended Darcy model still predicts slightly higher near-wall concentrations due to limited shear penetration, but the differences are smaller than at high $D_a$ . For $h = 0.5$ , this weak advection effect remains more visible than in a thinner layer (figure 15).

Increasing $h$ amplifies the impact of the Brinkman term at all Darcy numbers (figures 15, 16 and 17). For thick layers ( $h = 0.8$ ), shear penetration influences a large volume of the porous medium, producing pronounced differences between the models across the entire $D_a$ range. For moderate thickness ( $h = 0.5$ ), the differences are still substantial at high $D_a$ and remain visible at low $D_a$ . Thin layers ( $h = 0.3$ ) show the smallest – but still non-negligible – differences, particularly at high $D_a$ where $\ell _v \gtrsim h$ .

The extended Darcy (Brinkman) formulation is the physically consistent choice whenever (i) $D_a \gtrsim h^2$ or (ii) $h$ is large enough for shear penetration to influence dispersion. The classical Darcy model is only adequate for thin layers with low permeability, where shear effects are negligible.