Degradable Solute Transport in Porous Media with Variable Hydrodynamic Dispersion

1. Introduction

The analysis of solute transport in porous media is important for understanding and managing environmental processes such as groundwater pollution [1], nutrient leaching in agricultural soils [2,3,4], and the migration of pollutants in aquifers [5,6]. Porous media systems act as natural filters, but their effectiveness in controlling the dispersion of substances depends on the interaction between physical transport mechanisms and chemical-biological transformations and degradation [7]. In particular, the behavior of degraded solutes, including fertilizers, pesticides, and industrial pollutants, is of increasing interest due to their direct impact on soil and water quality, sustainable land management, and ecosystem conservation [8,9].

The main processes that affect the movement of solutes in porous media include hydrodynamic dispersion, convection, adsorption, and degradation [1,2,3,10,11]. Hydrodynamic dispersion, which combines molecular diffusion and mechanical mixing, governs the distribution of solute plumes under flow conditions. In contrast to idealized constant dispersion models, in natural subsurface environments the dispersion coefficient is variable and depends on the flow velocity, pore structure, and medium heterogeneity [12,13,14]. Taking into account the variability of dispersion is essential for adequate modeling of contaminant transport.

Furthermore, adsorption also plays a key role in regulating solute concentrations in porous media [11,15]. Through reversible or irreversible binding to the solid surface, adsorption can slow the transport of solutes and reduce peak concentrations. For degradable solutes, the interplay between adsorption and degradation is particularly important: adsorption can either slow the mobility of solutes, increasing degradation time, or protect pollutants from transformation processes, prolonging their existence in the environment.

In order to represent the combined effects of advection, dispersion, and degradation in porous media, researchers have recently developed sophisticated modeling techniques [15,16]. These models shed light on how pollutants in the environment proliferate, endure, or decompose in many natural systems. However, a lot of research still makes the oversimplified premise that dispersion is constant, which can either overestimate or underestimate the spreading and modification of solutes. This restriction presents difficulties for environmental management, especially when evaluating remediation initiatives, groundwater protection plans, and the long-term viability of soil-water systems.

By examining the transport of degradable solutes in porous media under circumstances of varying hydrodynamic dispersion, the current study fills this knowledge gap. This work helps to improve the depiction of contaminant migration in subterranean environments by creating and evaluating a mathematical model that takes into account both spatially variable dispersion and solute degradation. Policymakers, hydrogeologists, and environmental engineers involved in pollution remediation, sustainable agriculture, and groundwater protection will find value in the findings.

2. Problem Statement and Mathematical Model

A two-zone porous medium comprising active and passive regions is examined. The solute transport equation, considering adsorption and decay, can be expressed as [11,17]

$${\displaystyle \frac{\partial \left(\theta c\right)}{\partial t}}=-divJ-{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{\partial c_p}{\partial t}}-\theta \lambda c,$$

(1)

$$J=-\theta D\nabla c+\overrightarrow{w}c,$$

(2)

where c denotes the volumetric concentration, $c_a$ represents the concentration of the adsorbed substance in the active zone, $c_p$ signifies the concentration of the adsorbed substance in the passive zone, J indicates the flow density of the solute, $\theta$ refers to the porosity, $\lambda$ is the first-order decay coefficient (degradation), and $\overrightarrow{w}$ denotes the velocity. The terms ${\displaystyle \frac{\partial c_a}{\partial t}}$ and ${\displaystyle \frac{\partial c_p}{\partial t}}$ represent the adsorption processes in the active and passive zones, respectively, and pertain to the mass transfer between the liquid phase and the solid surface resulting from adsorption phenomena.

By substituting equation (2) into equation (1), one can see that:

$${\displaystyle \frac{\partial \left(\theta c\right)}{\partial t}}+div\left(\overrightarrow{w}c\right)=div(\theta D\nabla c)-{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{\partial c_p}{\partial t}}-\theta \lambda c.$$

(3)

The following can be derived from (3) in the case of a one-dimensional system

$${\displaystyle \frac{\partial \left(\theta c\right)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left(w_{x}c\right)={\displaystyle \frac{\partial }{\partial x}}\left(\theta D{\displaystyle \frac{\partial c}{\partial x}}\right)-{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{\partial c_p}{\partial t}}-\theta \lambda c,$$

(4)

where $w_x$ is the coordinate of the velocity along the $Ox$ axis.

If the medium is homogeneous, that is, $\theta =const$, the filtration rate is also constant, it is obtained that

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial x}}={\displaystyle \frac{\partial }{\partial x}}\left(D{\displaystyle \frac{\partial c}{\partial x}}\right)-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-\lambda c,$$

(5)

where $v_x={\displaystyle \frac{w_x}{\theta }}$ is the physical velocity of the fluid.

As mentioned above, in most cases hydrodynamic dispersion coefficient should be variable. As mentioned in [18] dispersion coefficient can be found in the form

$$D={\displaystyle \frac{1}{2}}{\displaystyle \frac{d{\sigma }^2}{dt}},$$

(6)

where $\sigma$ is the spatial variance of the tracer or solute distribution.

For longitudinal dispersivity ${\alpha }_L$

$${\alpha }_L={\displaystyle \frac{1}{2}}{\displaystyle \frac{d{\sigma }^2}{dx}},$$

(7)

$${\alpha }_L=f\left(x\right),$$

(8)

$${\sigma }^2=2\int f\left(x\right)dx.$$

(9)

Several suggested dispersivity function types (denoted linear, parabolic, asymptotic, and exponential) along with the corresponding theoretical variance functions obtained from (9) in [18] are

Linear:

$${\alpha }_L=ax,$$

(10)

$${\sigma }^2=a{x}^2.$$

(11)

Parabolic:

$${\alpha }_L=a{x}^b,$$

(12)

$${\sigma }^2=2a\left({\displaystyle \frac{x^{b+1}}{b+1}}\right).$$

(13)

Asymptotic:

$${\alpha }_L=A\left(1-{\displaystyle \frac{B}{x+B}}\right),$$

(14)

$${\sigma }^2=2A\left\{x-B\left[ln\left({\displaystyle \frac{x+B}{B}}\right)\right]\right\}.$$

(15)

Exponential:

$${\alpha }_L=E[1-exp(-Fx)],$$

(16)

$${\sigma }^2=2E\left[1+{\displaystyle \frac{exp\left(Fx\right)}{F}}\right],$$

(17)

where ${\alpha }_L$ longitudinal dispersivity, ${\sigma }^2$ spatial variance of solute distribution, x mean travel distance, $a,b,F$ constants, $A,E$ asymptotic or maximum dispersivity value, B characteristic half length (equals mean travel distance corresponding to $A/2$ ).

Hydrodynamic dispersion coefficient taken in following form in [12]

$$D=D_0+\epsilon \left|v\right|,$$

(18)

where $D_0$ represents the diffusion coefficient of a porous medium, which is frequently disregarded in field-scale tracer transport, and $\left|v\right|=\left|{\displaystyle \frac{q}{\theta }}\right|$ denotes the magnitude of the macroscopic pore water velocity and $\epsilon$ represent the dispersivity, which is traditionally regarded as a scale-invariant property of the porous medium.

These data indicate that longitudinal dispersivity, $\epsilon$, can be approximated by an empirical hyperbolic distance-dependent model of the following form

$$\epsilon ={\displaystyle \frac{1}{\frac{1}{{\epsilon }_{\infty }}+\frac{1}{{\beta }_x}}},$$

(19)

where ${\epsilon }_{\infty }$ is an asymptotic dispersivity [L] that is achieved at great distances, $\beta$ is a scale factor [L0] that describes the linear growth of the dispersion process when it is close to the origin, and x is the distance from the injection site.

In this paper used following expressions for hydrodynamic dispersion

$$D=D_0\left(1+e^{-\omega x}\right),D_0=const,\omega >0.$$

(20)

$$D=D_0\left(1+{\displaystyle \frac{x}{L}}\right),D_0=const.$$

(21)

$$D=D_0\left(1-{\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2}{L^2}}\right).$$

(22)

To determine the concentration of adsorbed substance in the active zone, following multistage kinetics [17,19] is used

$${\displaystyle \frac{\partial c_a}{\partial t}}=\left\{\begin{array}{cc}{\beta }_{ar}c-{\lambda }_a{c}_a,\hfill & 0<c_a\le c_{ar},\hfill \\ {\beta }_{aa}c-{\beta }_{ad}{c}_a-{\lambda }_a{c}_a,\hfill & 0<c_a<c_{ar},\hfill \\ 0,\hfill & c_a=c_{a0},\hfill \end{array}\right.$$

(23)

where $c_{a0}$ is the maximum concentration of adsorbed substance that can be achieved in the active zone, and ${\lambda }_a$ is the decay coefficient of the adsorbed substance in the liquid zone, The kinetic coefficients are represented by the variables ${\beta }_{ar}$, ${\beta }_{aa}$, and ${\beta }_{ad}$, whereas the variable $c_{ar}$ indicates the highest concentration at which the "charging" effect comes to a stop.

The following is the equation for the kinetics of adsorption in the passive zone [17,19]

$${\displaystyle \frac{\partial c_p}{\partial t}}=\left\{\begin{array}{cc}{\beta }_{p0}c-{\lambda }_p{c}_p,\hfill & 0<c_p\le c_{p1},\hfill \\ {\beta }_{p0}{\displaystyle \frac{c_{p1}}{c_p}}c-{\lambda }_p{c}_p,\hfill & c_{p1}<c_p<c_{p0},\hfill \\ 0,\hfill & c_p=c_{p0},\hfill \end{array}\right.$$

(24)

where ${\lambda }_p$ is coefficient of decay of the adsorbed substance in the passive zone, $c_{p1}$ is the concentration at which the "aging" effect begins.

Equations (5), (23), (24) are solved with the following initial and boundary conditions

$$c(0,x)=0,c_a(0,x)=0,c_p(0,x)=0.$$

(25)

$$c(t,0)=c_0=\mathrm{const},c(t,L)=0.$$

(26)

3. Numerical Solution

To solve problems (5), (23)-(26), with different expressions for hydrodynamic dispersion (20)-(22) we use the finite difference method [20,21]. In the domain $D=\{0\le x<\infty ,0\le t\le T\}$, we introduce a grid ${\omega }_{h\tau }$ , where T is the maximum time in the process under study. On the $OX$ axis, we divide the interval $[0,\infty )$ into pieces with a step of h and the interval $[0,T]$ into pieces with a step of $\tau$ along the time. To approximate the problem, we introduce the following grid

${\omega }_{h\tau }=\left\{\left(x_i,t_j\right),x_i=ih,i=0,1,\dots ,t_j=j\tau ,j=0,1,\dots ,J,\tau =T/J\right\}$.

Instead of the functions $c(t,x),c_a(t,x)$, $c_p(t,x)$, we consider the specific functions whose values at the nodes $\left(x_i,t_j\right)$ determine $c_i^j,c_{a,i}^j$, $c_{\mathrm{p},i}^j$, respectively.

In the case of for D used (20), in (5), we have

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial x}}={\displaystyle \frac{\partial }{\partial x}}\left[D_0\left(1+e^{-\omega x_i}\right){\displaystyle \frac{\partial c}{\partial x}}\right]-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c.$$

(27)

After some simplifications can be obtained

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial x}}=D_0\left(1+e^{-\omega x_i}\right)·{\displaystyle \frac{{\partial }^{2}c}{\partial x^2}}-D_0·\omega ·e^{-\omega x_i}·{\displaystyle \frac{\partial c}{\partial x}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c,$$

(28)

or

$${\displaystyle \frac{\partial c}{\partial t}}+\left(v_x+D_0·\omega ·e^{-\omega x_i}\right){\displaystyle \frac{\partial c}{\partial x}}=D_0\left(1+e^{-\omega x_i}\right)·{\displaystyle \frac{{\partial }^{2}c}{\partial x^2}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c.$$

(29)

Which can be approximated on the grid ${\omega }_{h\tau }$ in the following form

$$\begin{array}{c}{\displaystyle \frac{c_i^{j+1}-c_i^j}{\tau }}+\left(v_x+D_0·\omega ·e^{-\omega x_i}\right)·{\displaystyle \frac{c_i^{j+1}-c_{i-1}^{j+1}}{h}}=\hfill \\ \hfill D_0\left(1+e^{-\omega x_i}\right)·{\displaystyle \frac{c_{i+1}^{j+1}-2·c_i^{j+1}+c_{i-1}^{j+1}}{h^2}}-{\displaystyle \frac{1}{\theta }}·{\displaystyle \frac{c_{a,i}^{j+1}-c_{a,i}^{j+1}}{\tau }}-{\displaystyle \frac{1}{\theta }}·{\displaystyle \frac{c_{p,i}^{j+1}-c_{p,i}^{j+1}}{\tau }}-{\lambda }_e·c_i^{j+1}\end{array}$$

(30)

From which can be found

$$\begin{array}{c}{h}^2\left(c_i^{j+1}-c_i^j\right)+h·\tau ·\left(v_x+D_0·\omega ·e^{-\omega x_i}\right)·\left(c_i^{j+1}-c_{i-1}^{j+1}\right)=\tau ·D_0·\left(1+e^{-\omega x_i}\right)\hfill \\ \hfill \left(c_{i+1}^{j+1}-2·c_i^{j+1}+c_{i-1}^{j+1}\right)-{\displaystyle \frac{h^2}{\theta }}·{\displaystyle \frac{c_{a,i}^{j+1}-c_{a,i}^{j+1}}{\tau }}-{\displaystyle \frac{h^2}{\theta }}·{\displaystyle \frac{c_{p,i}^{j+1}-c_{p,i}^{j+1}}{\tau }}-{\lambda }_e·\tau ·h^2·c_i^{j+1}\end{array}$$

(31)

$$\begin{array}{c}\left(h·\tau ·\left(v_x+D_0·\omega ·e^{-\omega x_i}\right)+\tau ·D_0·\left(1+e^{-\omega x_i}\right)\right)·c_{i-1}^{j+1}-\left(h^2+h·\tau ·\left(v_x+\right.\right.\hfill \\ \left.+D_0·\omega ·e^{-\omega x_i}+2·\tau ·D_0\left(1+e^{-\omega x_i}\right)+{\lambda }_e·\tau ·h^2\right)·c_i^{j+1}+\tau ·D_0·\left(1+e^{-\omega x_i}\right)·c_{i+1}^{j+1}=\\ \hfill =-h^2·c_i^j+{\displaystyle \frac{h^2}{\theta }}·\left(c_{a,i}^{j+1}-c_{a,i}^j\right)+{\displaystyle \frac{h^2}{\theta }}·\left(c_{p,i}^{j+1}-c_{p,i}^j\right)\end{array}$$

(32)

Which can be written in following short form

$$A·c_{i-1}^{j+1}-B·c_i^{j+1}+E·c_{i+1}^{j+1}=-F_i^j,$$

(33)

where

$$A=h·\tau ·\left(v_x+D_0·\omega ·e^{-\omega x_i}\right)+\tau ·D_0·\left(1+e^{-\omega x_i}\right),$$

$$B=h^2+h·\tau ·\left(v_x+D_0·\omega ·e^{-\omega x_i}\right)+2·\tau ·D_0\left(1+e^{-\omega x_i}\right)+{\lambda }_e·\tau ·h^2,$$

$$E=\tau ·D_0·\left(1+e^{-\omega x_i}\right),$$

$$F=h^2·c_i^j-{\displaystyle \frac{h^2}{\theta }}·\left(c_{a,i}^{j+1}-c_{a,i}^j\right)-{\displaystyle \frac{h^2}{\theta }}·\left(c_{p,i}^{j+1}-c_{p,i}^j\right).$$

Equation (33) can be solved using Tridiagonal matrix algorithm [20,21]

In the case of for D used (21), in (5), we have

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left[D_0\left(1+{\displaystyle \frac{x}{L}}\right){\displaystyle \frac{\partial c}{\partial x}}\right]-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c,$$

(34)

or

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial t}}=D_0\left(1+{\displaystyle \frac{x}{L}}\right)·{\displaystyle \frac{{\partial }^{2}c}{\partial x^2}}-D_0·{\displaystyle \frac{1}{L}}·{\displaystyle \frac{\partial c}{\partial x}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c.$$

(35)

Which is approximated as

$$\begin{array}{c}{\displaystyle \frac{c_i^{j+1}-c_i^j}{\tau }}+\left(v_x+D_0{\displaystyle \frac{1}{L}}\right)·{\displaystyle \frac{c_i^{j+1}-c_{i-1}^{j+1}}{h}}=D_0\left(1+{\displaystyle \frac{x_i}{L}}\right)·{\displaystyle \frac{c_{i+1}^{j+1}-2·c_i^{j+1}+c_{i-1}^{j+1}}{h^2}}-\hfill \\ \hfill -{\displaystyle \frac{1}{\theta }}·{\displaystyle \frac{c_{a,i}^{j+1}-c_{a,i}^{j+1}}{\tau }}-{\displaystyle \frac{1}{\theta }}·{\displaystyle \frac{c_{p,i}^{j+1}-c_{p,i}^{j+1}}{\tau }}-{\lambda }_e·c_i^{j+1}\end{array}$$

(36)

Equation (36) can be solved as (30). Only difference will in equation (33) coefficients $A,B,E,F$ will have following form

$$A=h·\tau \left(v_x+D_0·{\displaystyle \frac{1}{L}}\right)+\tau ·D_0·\left(1+{\displaystyle \frac{x_i}{L}}\right),$$

$$B=h^2+h·\tau \left(v_x+D_0·{\displaystyle \frac{1}{L}}\right)+2·\tau ·D_0·\left(1+{\displaystyle \frac{x_i}{L}}\right)+{\lambda }_e·\tau ·h^2,$$

$$E=\tau ·D_0·\left(1+{\displaystyle \frac{x_i}{L}}\right),$$

$$F=h^2·c_i^j-{\displaystyle \frac{h^2}{\theta }}\left(c_{a,i}^{j+1}-c_{a,i}^j\right)-{\displaystyle \frac{h^2}{\theta }}\left(c_{p,i}^{j+1}-c_{p,i}^j\right).$$

In the case of for D used (22), in (5), we have

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left[D_0\left(1-{\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2}{L^2}}\right){\displaystyle \frac{\partial c}{\partial x}}\right]-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c,$$

(37)

or

$${\displaystyle \frac{\partial c}{\partial t}}+v_x{\displaystyle \frac{\partial c}{\partial t}}=D_0\left(1-{\displaystyle \frac{1}{2}}{\displaystyle \frac{x^2}{L^2}}\right)·{\displaystyle \frac{{\partial }^{2}c}{\partial x^2}}-D_0·{\displaystyle \frac{x}{L^2}}·{\displaystyle \frac{\partial c}{\partial x}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_a}{\partial t}}-{\displaystyle \frac{1}{\theta }}{\displaystyle \frac{\partial c_p}{\partial t}}-{\lambda }_{e}c.$$

(38)

Which is approximated as

$$\begin{array}{c}{\displaystyle \frac{c_i^{j+1}-c_i^j}{\tau }}+\left(v_x+D_0{\displaystyle \frac{x_i}{L^2}}\right)·{\displaystyle \frac{c_i^{j+1}-c_{i-1}^{j+1}}{h}}=D_0\left(1-{\displaystyle \frac{1}{2}}{\displaystyle \frac{x_i^2}{L^2}}\right)·{\displaystyle \frac{c_{i+1}^{j+1}-2·c_i^{j+1}+c_{i-1}^{j+1}}{h^2}}-\hfill \\ \hfill -{\displaystyle \frac{1}{\theta }}·{\displaystyle \frac{c_{a,i}^{j+1}-c_{a,i}^{j+1}}{\tau }}-{\displaystyle \frac{1}{\theta }}·{\displaystyle \frac{c_{p,i}^{j+1}-c_{p,i}^{j+1}}{\tau }}-{\lambda }_e·c_i^{j+1}\end{array}$$

(39)

In this case coefficients $A,B,E,F$ in equation (33) will have following form

$$A=h·\tau \left(v_x+D_0·{\displaystyle \frac{x_i}{L^2}}\right)+\tau ·D_0·\left(1-{\displaystyle \frac{x_i^2}{L^2}}\right),$$

$$B=h^2+h·\tau \left(v_x+D_0·{\displaystyle \frac{x_i}{L^2}}\right)+2·\tau ·D_0·\left(1-{\displaystyle \frac{1}{2}}{\displaystyle \frac{x_i^2}{L^2}}\right)+{\lambda }_e·\tau ·h^2,$$

$$E=\tau ·D_0·\left(1-{\displaystyle \frac{x_i^2}{L^2}}\right),$$

$$F=h^2·c_i^j-{\displaystyle \frac{h^2}{\theta }}\left(c_{a,i}^{j+1}-c_{a,i}^j\right)-{\displaystyle \frac{h^2}{\theta }}\left(c_{p,i}^{j+1}-c_{p,i}^j\right).$$

The differential schemes for equations (23) and (24) are as follows:

$${\displaystyle \frac{c_{a,i}^{j+1}-c_{a,i}^j}{\tau }}=\left\{\begin{array}{cc}{\beta }_{ar}\nu c_i^j-{\lambda }_{ea}{c}_{a,i}^j,\hfill & 0<c_{a,i}^j\le c_{ar},\hfill \\ {\beta }_{aa}v{c}_i^j-{\beta }_{a\mathrm{d}}{c}_{a,i}^j-{\lambda }_{ea}{c}_{a,i}^j,\hfill & c_{ar}<c_{a,i}^j<c_{a0},\hfill \\ 0,\hfill & c_{a,i}^j=c_{a0}.\hfill \end{array}\right.$$

(40)

$${\displaystyle \frac{c_{p,i}^{j+1}-c_{p,i}^j}{\tau }}=\left\{\begin{array}{cc}{\beta }_{p0}{c}_i^j-{\lambda }_{ep}{c}_{p,i}^j,\hfill & 0<c_{p,i}^j\le c_{p1},\hfill \\ {\displaystyle \frac{{\beta }_{p0}{c}_{p1}{c}_i^j}{c_{p,i}^j}}-{\lambda }_{ep}{c}_{p,i}^j\hfill & c_{p1}<c_{p,i}^j<c_{p0},\hfill \\ 0,\hfill & c_{p,i}^j=c_{p0}.\hfill \end{array}\right.$$

(41)

Initial and boundary conditions (25) and (26) have the form

$$c_{a,i}^j=0,i=\overline{0,I},j=0,c_{p,i}^j=0,i=\overline{0,I},j=0,$$

(42)

$$c_i^j=0,i=\overline{0,I},j=0,c_i^j=c_0,i=0,j=\overline{0,J}.$$

(43)

4. Results and Discussion

To carry out numerical experiments we developed a program in Python.

Figure 1 shows the results for the case where $D_0=0$, i.e., diffusion and hydrodynamic dispersion are not taken into account. The results are presented in the form of concentration profiles $c/c_0,c_a$ and $c_p$. Over time, the values of $c/c_0,c_a$ and $c_p$ at fixed points in the reservoir increase (Figure 1). It can be seen from the results that the solute transport process occurs more slowly when diffusion is not taken into account. Figure 1.a shows that at $t=3000\mathrm{s}$ the substance concentration has only spread over a distance of $x\approx 0.32\mathrm{m}$. At this time, the adsorption in both zones has approached $x\approx 0.3\mathrm{m}$.

Figure 2 shows the case where $D_0=5·{10}^{-6}$ is constant diffusion. Comparing with Figure 1, it can be concluded that taking into account diffusion significantly accelerates the solute transport and adsorption. In particular, at $t=3000\mathrm{s}$ in Figure 2.a, it can be observed that the substance concentration reaches the end of the medium ( $L=0.4\mathrm{m}$ ). This, in turn, is observed with an increase in the concentration of the adsorbed substance in both zones (Figure 2b, Figure 2c).

Figure 3 presents comparative graphs for different values of $D_0$. It can be seen from the graphs that an increase in the value of diffusion leads to a wider spread of the concentration profiles towards the interior of the medium. Since the scattering is less at $D_0=0$, the concentration $c/c_0$ at points closer to the point $x=0$ is greater than at the point $D_0\ne 0$, and vice versa at points further from the point $x=0$. This was not observed for the concentrations $c_a$ and $c_p$, i.e. their values increased significantly at all points of the medium with increasing diffusion.

Figure 4 presents the results of numerical experiments conducted according to formula (20), where an exponential expression depending on the coordinate for the dispersion coefficient is obtained. The results show that using formula (20) instead of the expression $D=D_0$ increases the concentration distribution, albeit partially. This, in turn, increases the amount of adsorbed substance in both the active and passive zones.

In Figure 5 shown the solute transport and adsorption changes for different values of $\omega$ ($\omega =0,\omega =1,\omega =10$). It can be seen from the graphs that the cases with $\omega =0$ and $\omega =1$ do not differ much from each other. When $\omega =0$, it can be observed that all three concentrations increase by a very small amount compared to when $\omega =1$. However, when $\omega =10$, it can be observed that the concentrations decrease sharply at points further from the point $x=0$.

Figure 6 illustrates the outcomes of numerical experiments performed in accordance with formula (21), resulting in a linear increasing expression dependent on the coordinate for the dispersion coefficient. The results indicate that employing formula (21) in place of the expression $D=D_0$ leads to a minimal increase in the concentration distribution, which may not be immediately apparent upon initial observation.

Figure 7 presents the results of numerical experiments conducted in alignment with formula (22), yielding a parabolic decreasing function related to the coordinate for the dispersion coefficient. The findings demonstrate that utilizing formula (22) instead of the expression $D=D_0$ yields results that are highly comparable.

Figure 8 compares the results obtained for expressions (20), (21), (21) at $t=9000\mathrm{s}$. It can be seen from the results that the results obtained for expressions (21) and (22) do not differ much. In the results obtained for expression (20), the values of $c_a$ and $c_p$ are significantly larger. The values of $c/c_0$ are smaller at points closer to the point $x=0$ for expression (20) than for (21) and (22), and vice versa, the distance at $x=0$ is larger.

Comparative Table 1, Table 2 and Table 3 are presented for a more in-depth analysis of the results obtained using the different expressions proposed for dispersion. It can be seen from Table 1 that when the expression (20) is used, the concentrations at all fixed points are larger in the cases of $\omega =0$ and $\omega =1$ than when $D=D_0$ is constant. When $\omega =10$, the difference is not so great. When the expression (21) is used, there is a difference, although it is small compared to $D=D_0$, but for the expression (22) the difference is almost imperceptible in the results obtained with 4 digits of accuracy. This can be explained by the fact that due to the small values of x, its square becomes even smaller.

Tables Table 2, Table 3 present the results for larger values of time, and the conclusions made for Table 1 are once again confirmed here.

5. Conclusions

This study formulated and examined a degradable-solute transport model inside a two-zone porous medium that explicitly integrates advection, spatially-dependent hydrodynamic dispersion, multistage adsorption in active and passive zones, and first-order decay. A finite-difference solution method utilizing a tridiagonal solver was employed to address coordinate-dependent diffusion terms and piecewise adsorption kinetics. Numerical tests were conducted for constant dispersion and three representative spatially variable forms—exponential, linear, and parabolic—using identical hydraulic and kinetic parameters.

The simulations indicate that disregarding diffusion and dispersion significantly hinders plume progression and leads to an underestimation of both aqueous concentrations and sorbed inventory. The introduction of even minimal dispersion enhances transport throughout the domain and amplifies adsorption in both areas, resulting in reduced aqueous peaks at the input and elevated concentrations further downstream. Among the variable-dispersion forms, the exponential law dispersion exhibited the most significant deviations from the constant dispersion scenario. Conversely, the linear and parabolic dispersion forms produced minimal variations at the examined length and velocity scales, indicating their relatively modest variation within specified boundaries. The sensitivity to the exponential-form parameter $\omega$ indicated that a quick decay of the elevated inlet dispersion ($\omega =10$) inhibits far-field accumulation while just slightly influencing the near-inlet profile. These results collectively indicate that the shape and scale of the dispersion function govern both the timing of breakthrough and the distribution between dissolved and sorbed masses.

The findings suggest that utilizing scale-dependent or coordinate-dependent dispersion models can significantly impact forecasts of pollutant arrival timings and retention in remedial or agricultural contexts. Augmented dispersion near sources typically results in a more extensive distribution of solute, hence enhancing sorptive absorption and potentially reducing local peak water concentrations; nevertheless, it may also expedite downstream exposure if the greater dispersion endures across significant distances. Consequently, depending on a fixed dispersion coefficient may inaccurately represent danger and cleanup timelines in systems where dispersion fluctuates with distance or hydraulic conditions.

This study is confined to one-dimensional, steady-flow scenarios with uniform porosity and simplified boundary conditions, excluding calibration to laboratory or field data. Future endeavors should broaden the framework to encompass two and three dimensions with heterogeneous characteristics, integrate it with variable-density or transient flows, investigate alternative (nonlinear) adsorption/desorption and biodegradation kinetics, and conduct systematic sensitivity and uncertainty evaluations. Integrating data-driven calibration or inverse modeling with tracer experiments would enhance the quantification of the suitable form and parameters of $D\left(x\right)$ for site-specific forecasts. Notwithstanding these constraints, the study unequivocally demonstrates that the careful selection of dispersion models is crucial for accurate predictions of degradable-solute behavior and movement in porous surfaces.