A Peridynamic Differential Operator-Based Model for Quantifying Spatial Non-Local Transport Behavior of Pollutants in Heterogeneous Media

Abstract

Modeling pollutant transport in heterogeneous media is an important task of hydrology. Pollutant transport in a non-homogeneous environment typically exhibits non-local transport dynamics, whose efficient characterization requires a parsimonious model with the non-local feature. This study encapsulates the non-local transport characteristic of pollutants into the peridynamic differential operator (PDDO) and develops a PDDO-based model for quantifying the observed pollutant non-local transport behavior. The simulation results show that the proposed model can describe pollutant non-local transport behavior in various heterogeneous media. The non-local nature of pollutant transport can be adjusted by pre-defined weight function $w\left(\right|\xi \left|\right)$ and horizon $H_x$. Applications show that the PDDO-based model can better capture pollutant non-local transport behavior than the classical advection–diffusion equation (ADE) model, especially for quantifying the tail of the experimental data late. Analyses further reveal that the PDDO-based model can characterize both normal (Fickian) and anomalous (Lévy) diffusion regimes.

1. Introduction

The efficient modeling of pollutant transport in a natural environment (surface and subsurface water) is one of the major issues of computational hydrology [1,2,3,4,5]. Pollutant transport in heterogeneous media exhibits complex dynamics, due to the multi-scale heterogeneity of the media, the complex exchange processes between water and tracer, and complex turbulence in the natural environment [5,6,7]. A key problem is how to accurately characterize pollutant transport with complex dynamics in heterogeneous media. Thus, many models which are capable of well describing the complex processes of pollutant transport in heterogeneous media have attracted extensive attention in the last two decades, especially the non-local models [8,9,10,11]. It is also noteworthy that although the physical heterogeneity of the porous medium is essential, many studies have focused on the effects of spatially variable geochemical heterogeneities because the heterogeneous soil cannot be treated as uniform and homogeneous in the natural subsurface environment. The spatially variable geochemical heterogeneities affect pollutant transport significantly [12,13,14,15]. Additionally, many researchers have investigated the transport behavior of molecular size pollutants as well as suspended particle pollutants and found that the transport behavior between the solute type and suspended particle type of pollutants is very different [15,16]. In this study, we focus on the non-local pollutant transport in heterogeneous media, which is usually characterized by an apparent early arrival feature of the pollutant plume and can pose a high risk to the environment [17]. The dynamics of pollutant non-local transport in the heterogeneous systems are not fully understood, motivating this study.

The non-local transport dynamics of pollutant particles may be related to the preferential flow paths in heterogeneous systems [18]. As shown in Figure 1, the pollutant particles may experience a sudden entrainment process, and they can travel considerably far along the “preferential flow paths”. Considering the non-local transport behaviors of pollutants, various non-local models have been proposed [19,20,21]. The main discrepancy between the non-local and classical local models is that non-local models usually refer to the integro-differential expression to characterize the spatial interaction natures in heterogeneous systems. The non-local transport behavior of pollutants exists at all scales: from the macroscale in field experiments to laboratory experiments, especially on the micro-scale. Hence, the existing models in prediction are sometimes acceptable in practice, but not always. The measured concentration flux can differ by one to two orders of magnitude from predictions [22].

However, despite considerable efforts over the last two decades, the capability to accurately characterize pollutant transport behavior in the heterogeneous system remains inadequate. One of the major obstacles is the poor understanding and modeling of non-local transport behaviors of pollutants in a wide range of conditions. The advection–dispersion equation (ADE) model provides a classical framework for characterizing pollutant transport and makes an important contribution to the following theoretical development. However, a growing number of researchers have found that the pollutant moving through non-homogeneous media (surface or subsurface water) cannot be described by using the classical ADE model. This study aims to quantify pollutant non-local transport behaviors in a wide range of conditions and proposes a peridynamic differential operator (PDDO) based model for quantifying the observed non-local transport behavior. PDDO was first built by Madenci et al. [23,24], who proposed the PDDO by introducing the classical peridynamic (PD) theory [25]. The PDDO can describe any order of partial derivatives of the spatial and temporal functions by using the orthogonality property of the PD functions. The PDDO is expressed in terms of only integration over the domain of the pre-defined horizon. Hence, PDDO converts the local differentiation equation to its non-local integration form in a unified manner, and it is not prone to the singularities that arise from the presence of discontinuities in modeling [24,26]. Here, this study introduces a PDDO into the ADE model to describe the dynamics of pollutant non-local transport. Considering various pre-defined weight functions and horizon domains, the proposed PDDO-based model can characterize pollutant non-local transport behavior in various conditions. The classical local model is extended to a general region scale model, even on a global scale, as shown in Figure 2.

The literature reviewed above takes into account the basic concept maps of the pollutant non-local transport and lays the foundation for the wide application of the PDDO-based model. The purpose of this study is to answer the following three questions using the proposed PDDO-based model. First, can the PDDO-based model bridge the pollutant transport’s local and regional scales? Second, how do the pre-defined weight function and horizon domains of the PDDO-based model affect the dynamics of the pollutant non-local transport process? Third, can the PDDO-based model accurately characterize the spatial correlation of pollutant non-local transport?

The rest of this work is organized as follows. Section 2 reviews the PDDO theory and proposes a PDDO-based model. Section 3 gives an analytical solution and numerical algorithm of the PDDO-based model. One-dimensional and two-dimensional results of the PDDO-based model with different weight functions and horizon domains are calculated, and the applications of the proposed non-local model are also shown in this section. Section 4 discusses the diffusion regime of the PDDO-based model with various PD functions, and the limitations of the PDDO-based model are also stated here. Conclusions are finally drawn in Section 5.

2. Model Development

Pollutant transport simulation is one of the core tasks of computational hydrology, and various physical models have been developed and applied in the past decades [21,27,28,29]. Mathematical modeling and understanding pollutant transport requires determining spatial derivatives of the concentration. However, many studies have found that pollutants exhibit non-local transport behavior, especially in non-homogeneous environments, which cannot be well characterized by the classical local definition of spatial derivatives in modeling [17,30]. Considering the drawbacks of the differentiation processes, the integration processes are more suitable for characterizing non-local pollutant transport. Integration is a typical non-local process because it depends on the entire range of integration. Hence, the PDDO is introduced in this study to characterize the pollutant non-local transport behavior.

The most widely used model ADE is always applied to characterize the pollutant transport behavior, which, however, cannot well describe the non-local transport behavior of pollutants in heterogeneous media [21,31]. The classical ADE types model can be expressed as follows [32,33,34,35,36]:

$$\begin{array}{c}\hfill R\frac{\partial c(x,t)}{\partial t}=-v\frac{\partial c(x,t)}{\partial x}+D\frac{{\partial }^{2}c(x,t)}{\partial x^2}\end{array}$$

(1)

where $c(x,t)$ is a representation of the pollutant concentration as a function of the space and time, v [L/T] and D [L${}^2$/T] are the mean advection velocity and diffusion coefficient of pollutant particles, respectively. R (dimensionless) is the retarded coefficient. Thus, Equation (1) is called the retarded advection–diffusion equation (R-ADE) [37].

The above ADE-type models have been widely used in many studies to describe the pollutant transport process. However, the observed non-local pollutant transport cannot be precisely described by the R-ADE model. The non-local transport that exists in various systems proposes, as do many studies, the following fractional ADE (F-ADE) model to characterize the observed non-local transport behavior of pollutants:

$$\begin{array}{c}\hfill R\frac{\partial c(x,t)}{\partial t}=-v\frac{\partial c(x,t)}{\partial x}+D\frac{{\partial }^{\beta }c(x,t)}{\partial x^{\beta }}\end{array}$$

(2)

where ${\partial }^{\beta }/\partial x^{\beta }$ is the fractional operator. In this study, the Riemann–Liouville fractional derivative is considered in the simulation [38].

The F-ADE model provides a robust tool to describe the non-local pollutant transport along the preferential flow path, as shown in Figure 1. The basic form of space F-ADE replaces the integer-order derivative in space with a fractional one $\beta$ (1<$\beta \le$2). When $\beta =2$, the F-ADE model reduces to the classic ADE. However, various environmental factors can affect pollutant transport in natural conditions. The F-ADE model cannot capture the non-locality of the bedload transport at various spatial scales due to the mathematical definition of the fractional operator. This theoretical deficiency limits researchers from further understanding the pollutant non-local transport behavior. Large deviations can arise from the lack of comprehensive quantification of the pollutant non-local transport process. Hence, we will propose a non-local model for quantifying pollutant non-local transport behavior with various conditions.

2.1. A Brief Review of PDDO

Unlike the classical methods (such as finite difference method, finite element method, etc.) for solving the advection–diffusion equation (ADE) type models to characterize the pollutant transport, the PDDO can be used to reconstruct the integral form of the target differential equation [24,26]. Hence, PDDO is an alternative tool to describe the solute non-local transport. Although the differentiation process is usually more direct than integration in analytical mathematics, the reverse is true in computational mathematics, especially in the presence of a jump discontinuity or a singularity. Integration is a non-local process because it depends on the entire range of integration. Considering non-local pollutant transport, the PDDO is an efficient tool for describing the pollutant non-local transport behavior along the preferential paths.

The PDDO can be obtained by using the Taylor series expansion of a function $f(x+\xi )$ as follows:

$$\begin{array}{c}\hfill f(x+\xi )-f\left(x\right)=\xi ·\frac{df\left(x\right)}{dx}+\frac{{\xi }^2}{2}·\frac{d^{2}f\left(x\right)}{d{x}^2}+\dots +\frac{{\xi }^N}{N!}·\frac{d^N}{d{x}^N}f\left(x\right)+R\left(x\right)\end{array}$$

(3)

where $R\left(x\right)$ represents the remainder.

Here, a orthogonality property of PD functions, $g_N^p\left(\xi \right)$, is introduced, which satisfies the following equations [23]:

$$\begin{array}{c}\hfill \frac{1}{n!}\underset{H_x}{\int }{\xi }^n{g}_N^p\left(\xi \right)d\xi ={\delta }_{np}\mathrm{with} (n,p=0,1,\dots ,N)\end{array}$$

(4)

and

$$\begin{array}{c}\hfill g_N^p\left(\xi \right)=\sum _{q=0}^N{a}_q^p{w}_q\left(\xi \right){\xi }^q\end{array}$$

(5)

where $w_q\left(\xi \right)$ is the weight functions, and it can be used to measure the strength of the influence of the integral points in the defined horizon $H_x$.

Assuming that the value of the remainder, $R\left(x\right)$, is negligibly small, and multiplying each term in Equation (3) by the proposed PD functions $g_N^p\left(\xi \right)$ and integrating over the family of point, x, defined as $H_x=x^{\prime }\in [a,b]$, one obtains

$$\begin{array}{c}\hfill \underset{H_x}{\int }\left[f(x+\xi )-f\left(x\right)\right]d\xi =\frac{\partial f\left(x\right)}{\partial x}\underset{H_x}{\int }\xi g_N^p\left(\xi \right)d\xi +\frac{{\partial }^{2}f\left(x\right)}{\partial x^2}\underset{H_x}{\int }\frac{1}{2}{\xi }^2{g}_N^p\left(\xi \right)d\xi +\dots \\ \hfill +\frac{{\partial }^{N}f\left(x\right)}{\partial x^N}\underset{H_x}{\int }\frac{1}{N!}{\xi }^N{g}_N^p\left(\xi \right)d\xi +R(N,x)\end{array}$$

(6)

It is worth noting that each point has its family members in the range of the horizon $H_x$ and occupies an infinitesimally small entity. Hence, each family’s properties (such as size, shape, etc.) can be very different, and they significantly affect the strength of non-locality. The strength of the interaction between the family points in each family group is specified by the proposed weight functions $w_q\left(\xi \right)$.

Based on the orthogonal nature of the function $g_N^p\left(\xi \right)$, the spatial derivatives can be expressed as follows:

$$\begin{array}{c}\hfill \frac{d^{p}f\left(x\right)}{d{x}^p}=\underset{H_x}{\int }\left[f(x+\xi )-f\left(x\right)\right]g_N^p\left(\xi \right)d\xi \end{array}$$

(7)

Equation (7) gives a integration expressions of the spatial derivatives in the horizon $H_x$. In the regular PDDO, the shape of the horizon is specified as a sphere (3 dimensions), cylinder (2 dimensions), and segment (1 dimension) of the modeling. Further, the characteristic internal length parameter (radius of the defined sphere, circle, or line), $\delta$, referred to as the “horizon”, is constant.

2.2. PDDO-Based Model

As described above, pollutants may exhibit non-local transport behavior in heterogeneous media, and the classical ADE types model (1) cannot well characterize the non-local transport behavior. Considering the drawbacks of the model (1), the PDDO is introduced here to characterize the non-local transport behavior of pollutants. Figure 2 plots the conceptual map of the defined horizon $H_x$, and the local field of the model (1) is drawn for comparison. As shown in Figure 2, the horizon $H_x$ determines the non-local transport zone of pollutants, which also indicates the length of the preferential path in the media. A longer preferential path implies a stronger non-local transport behavior of pollutant particles, resulting in a larger horizon $H_x$.

By introducing the PDDO into Equation (1), the following PDDO-based non-local transport model can be obtained:

$$\begin{array}{c}\hfill R\frac{\partial c(x,t)}{\partial t}=\underset{H_x}{\int }\left(-v·g_N^1\left(\xi \right)+D·g_N^2\left(\xi \right)\right)\left[c(x+\xi ,t)-c\left(x\right)\right]d\xi \end{array}$$

(8)

the parameters v, D, and R are same as Equation (1), and the PD functions $g_N^1\left(\xi \right)$ and $g_N^2\left(\xi \right)$ are determined by using Equations (4) and (5).

Since $g_N^p\left(\xi \right)$ is zero outside of the defined horizon $H_x$, Equation (8) can be written as follows:

$$\begin{array}{c}\hfill R\frac{\partial c(x,t)}{\partial x}=-v·\left[g_N^1\ast c\right](x,t)+D·\left[g_N^2\ast c\right](x,t)-\beta c(x,t)\end{array}$$

(9)

where $(\ast )$ denotes the convolution operator:

$$\begin{array}{c}\hfill \left[g_N^p\ast c\right](x,t)={\int }_R{g}_N^p\left(\xi \right)c(x+\xi ,t)d\xi \end{array}$$

(10)

and

$$\begin{array}{c}\hfill \beta =\underset{H_x}{\int }\left(-v·g_N^1\left(\xi \right)+D·g_N^2\left(\xi \right)\right)d\xi \end{array}$$

(11)

Thus, Equation (9) is the PDDO-based model for characterizing the non-local transport behavior of pollutants. The analytical solution and numerical algorithm for solving the PDDO-based model (9) will be further given in the following section.

2.3. Numerical Algorithm of the PDDO-Based Model

As described above, the analytical form of the PDDO-based model cannot be obtained in most cases. Hence, the numerical algorithm is herein developed to solve the proposed PDDO-based model.

The one-dimensional case of the PD functions $g_N^p\left(\xi \right)$ can be set as follows [39]:

$$\begin{array}{c}\hfill g_2^p\left(\xi \right)=a_1^p{w}_1\left(\left|\xi \right|\right)\xi +a_2^p{w}_2\left(\left|\xi \right|\right){\xi }^2,p=1,2\end{array}$$

(12)

Considering the same weigh functions $w_1\left(\left|\xi \right|\right)$ and $w_2\left(\left|\xi \right|\right)$, for all terms, the Equation (4) can be expressed as follows:

$$\begin{array}{c}\hfill \sum _{q=1}^2{A}_{nq}{a}_q^p=b_n^p\end{array}$$

(13)

where

$$\begin{array}{c}\hfill A_{nq}=\underset{H_x}{\int }w\left(\left|\xi \right|\right){\xi }^{n+q}d\xi \end{array}$$

(14)

Equation (13) can be written as follows:

$$\begin{array}{c}\hfill \mathbf{Aa}=\mathbf{b}\end{array}$$

(15)

and

$$\begin{array}{c}\hfill \mathbf{a}={\mathbf{A}}^{-1}\mathbf{b}\end{array}$$

(16)

where

$$\begin{array}{c}\hfill \mathbf{A}=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\end{array}$$

(17)

and

$$\begin{array}{c}\hfill \mathbf{a}=\left[\begin{array}{cc}{a}_1^1& a_1^2\\ a_2^1& a_2^2\end{array}\right],\mathbf{b}=\left[\begin{array}{cc}1& 0\\ 0& 2\end{array}\right]\end{array}$$

(18)

Inserting $\mathbf{a}$ into Equation (12), $g_2^p\left(\xi \right)$ can be obtained. The implicit discrete form of Equation (8) can be expressed as follows:

$$\begin{array}{c}\hfill R\frac{c_i^{j+1}-c_i^j}{\Delta t}=\sum _{k=1}^{N_i}\left(-v·g_2^1\left(\xi \right)+D·g_2^2\left(\xi \right)\right)\left[c_k^{j+1}-c_i^{j+1}\right]\Delta x\end{array}$$

(19)

where $N_i$ is the number of the family point of the i-th space point, $\Delta x$ and $\Delta t$ are space interval and time interval, respectively.

Equation (19) gives an implicit form of the proposed PDDO-model. Compared with the discrete form of the classical ADE models by using the finite difference method, the discrete form of the PDDO-based model also indicates that the target point is not only affected by the nearby points, but also affected by all family points in the defined near-field $H_x$.

3. Results

3.1. Analytical Solution of PDDO Model

PDDO-based model (9) gives a powerful tool to solve the non-local transport behavior of pollutants, and the PD functions $g_N^p\left(\xi \right)$ and horizon $H_x$ determine the strength of the non-local property of pollutant transport. The horizon $H_x$ and weight function $w\left(\xi \right)$ are pre-defined according to the environment state before the simulation. In the simulation, the PD functions $g_N^p\left(\xi \right)$ should be analytically or numerically determined according to the pre-defined horizon $H_x$ and weight function $w\left(\xi \right)$, which are the core tasks in modeling.

In order to analytically solve the proposed PDDO-based model (9), taking the Fourier to transform Equation (9), one obtains

$$\begin{array}{c}\hfill R\frac{\partial \widehat{c}(k,t)}{\partial t}=-v·\left({\widehat{g}}_N^1\left(k\right)·\widehat{c}(k,t)\right)+D·\left({\widehat{g}}_N^2\left(k\right)·\widehat{c}(k,t)\right)-\beta \widehat{c}(k,t)\end{array}$$

(20)

where k is the Fourier variable.

Collating Equation (20) yields

$$\begin{array}{c}\hfill \frac{\partial \widehat{c}(k,t)}{\partial t}=\left\{\frac{1}{R}\left[-v·{\widehat{g}}_N^1\left(k\right)+D·{\widehat{g}}_N^2\left(k\right)-\beta \right]\right\}\widehat{c}(k,t)\end{array}$$

(21)

Considering the initial point source condition $\widehat{c}(k,0)=1$, the analytical solution of Equation (21) is

$$\begin{array}{c}\hfill \widehat{c}(k,t)=exp\left(\frac{1}{R}\left[-v·{\widehat{g}}_N^1\left(k\right)+D·{\widehat{g}}_N^2\left(k\right)-\beta \right]·t\right)\end{array}$$

(22)

Taking inverse Fourier transform of Equation (22) yields

$$\begin{array}{c}\hfill c(x,t)=F^{-1}\left\{exp\left(\frac{1}{R}\left[-v·{\widehat{g}}_N^1\left(k\right)+D·{\widehat{g}}_N^2\left(k\right)-\beta \right]·t\right)\right\}\left(k\right)\end{array}$$

(23)

Equation (23) gives a general form of the analytical solution of the proposed PDDO-based model (9). In order to obtain an explicit expression of the analytic solution (23), it is necessary to obtain the form of the function $g_2^p\left(\xi \right)$ in the Fourier space. However, it is impossible to obtain an analytical form of the function $g_2^p\left(\xi \right)$ in the Fourier space in most cases. Therefore, we will use the numerical algorithm to solve the PDDO-based model (9), as described in the section methodology.

3.2. One-Dimensional Results of the PDDO-Based Model

PDDO-based model gives an alternative tool to describe the non-local transport behavior of pollutants, and the non-local property of pollutants are determined by the weight function $w\left(\left|\xi \right|\right)$ and $w_2\left(\left|\xi \right|\right)$ and horizon $H_x$. To investigate the effect of the weight function on the pollutant transport process, the following power-law form of the weight function is considered:

$$\begin{array}{c}\hfill w\left(\left|\xi \right|\right)=\frac{{\delta }^k}{{\left|\xi \right|}^k}(k=0,1,2,3\dots )\end{array}$$

(24)

Figure 3 shows the weight function with different parameter k, and an exponential function $exp(-{\left|\xi \right|}^2/{\delta }^6)$ is also drawn for comparison. The result shows that the strength of point influence in the horizon increases as k decreases, and a smaller k implies a stronger non-local nature of pollutant transport. The exponential form of the weight function exhibits the weakest non-local property to characterize pollutant transport.

Figure 4 plots the breakthrough curves (BTCs) and the snapshots of the PDDO-based model (9) with different weight functions; the results of the classical ADE model are also drawn for comparison. As shown in Figure 4a, the PDDO-based model gives a significant leading edge characteristic in the early time, which means that the pollutants described by the non-local model can reach the observation site much earlier than described by the local model. Furthermore, the early arrival characteristic of pollutants becomes more evident with the enhancement of the non-local nature, and this is mainly because the pollutants are easier to reach the observation point through the preference paths. Figure 4b shows the snapshots of the PDDO-based model, and the results indicate that the non-local transport behavior of pollutants exhibits trailing edge characteristics. Compared with the ADE model, the snapshots of the PDDO-based model show a significantly asymmetric characteristic. The difference between the PDDO model and the ADE model is that the non-local nature is considered in the PDDO model. The concentration of the pollutant is not only locally relevant but also correlates with regions within the non-local area. Therefore, the introduction of the non-local nature leads to pollutant particles with early arrival features in BTCs and a significant trailing edge in the spatial snapshots.

3.3. Two-Dimensional Results of the PDDO-Based Model

Equation (9) gives a non-local description of one-dimensional pollutant transport. To further investigate the non-local transport behavior of pollutants in heterogeneous media, here we give the following two-dimensional form of the PDDO-based model. We assume that the advection process occurs only in the x-direction. One obtains

$$\begin{array}{c}\hfill R\frac{\partial c(x,y,t)}{\partial t}=\underset{H_x}{\int }\left[-v·g_2^{10}\left(\xi \right)+D_x·g_2^{20}\left(\xi \right)+D_y·g_2^{02}\left(\xi \right)\right]\\ \hfill ·\left(c(x+{\xi }_1,y+{\xi }_2,t)-c(x,y,t)\right)d{V}_x\end{array}$$

(25)

where $g_2^{10}$, $g_2^{20}$, and $g_2^{02}$ are PD functions; see Madenci et al. [24] for the detailed derivation.

The last subsection discusses the weight function for characterizing the pollutant non-local transport. Here, we compared two different horizon ranges $H_x$ in describing pollutant transport. Figure 5 shows the simulation results of the two-dimensional PDDO-based model (25); two different $H_x$ are considered, i.e., $H_x=0.6$$\Omega$ and $H_x=0.05$$\Omega$. The results indicate that the plume gives a symmetric (Gaussian) shape in a smaller $H_x$. This is mainly because the PDDO-based model can regress to a local ADE form, and the analytical solution of the ADE is Gaussian distribution in space.

It is worth noting that the domain of the PDDO-based model is divided into a finite number of cells, and the interior points can be assigned a symmetric family in space, while the points near the boundary have non-symmetric families. Therefore, the points near the boundary have their own PD function $g_N^p\left(\xi \right)$. As shown in Figure 5a, the large horizon $H_x$ give a large number of non-symmetric families of the points, and the plumes show significantly asymmetric features.

3.4. Applications

Here, we check the applicability of the PDDO-based model (9) in quantifying pollutant non-local transport documented in the literature [17]. To explore the dynamics of tracers transport in the alluvial setting with heterogeneity, Yin et al. [17] built a two-dimensional, different alluvial setting with various hydrofacies structures generated by using the T-PROGS, and simulated the conservative tracer transport by using the Monte-Carlo method. The detailed information can be found in the literature [17].

Hence, Yin et al. [17] gives an ideal set of data for investigating the non-local pollutant transport in heterogeneous media. As shown in Figure 6, four snapshots of the tracers are documented in the studies, i.e., 27 days, 132 days, 224 days, and 328 days. The experimental data show that the trailing edge becomes stronger with time evolution, which indicates that the non-local transport behavior enhances with time evolution.

We used the first snapshot (27 days) to determine the parameters v, D, and R for the PDDO-based and F-ADE models and then predicted the later snapshots by adjusting the parameters $H_x$ and $\beta$. Considering the parameters v and D of the R-ADE model have the same dimensions as those of the PDDO-based model, the parameters v and D in the R-ADE model were set to be the same as those of the PDDO-based model, andonly the parameter R was adjusted. The results of the ADE and F-ADE models are calculated for comparison. The weight function of the PDDO are pre-defined as $w\left(\right|\xi \left|\right)={\delta }^2/{\left|\xi \right|}^2$ in the simulations. Figure 6 plots the fitting results of the PDDO-based model for each snapshot. The results show that the PDDO-based and F-ADE models can well capture the non-local transport behavior, while the ADE model cannot capture the observation results. Moreover, the F-ADE model can capture the trailing edge of the snapshot well, but it is not accurate for the peak position. This phenomenon may be due to the F-ADE model being a global non-local model, which assumes that the pollutant transport process is affected by the entire space, leading it to overestimate the non-local transport behavior of the pollutant. The main discrepancy between the PDDO-based model and the ADE model is that the PDDO-based model can well capture the trailing edge feature of the experimental data. The ADE model is a local definition, which cannot capture the observed non-local pollutant transport behavior. Additionally, the experimental data show significant asymmetric features in the snapshots, and the ADE can describe the symmetric transport behavior of the pollutants.

Table 1. The best-fit parameters for F-ADE ($v_1$, $D_1$, $R_1$, $\beta$), ADE ($v_2$, $D_2$, $R_2$) and PDDO-based ($v_2$, $D_2$, $R_3$, $H_x$) models.

Time (Days)	${\mathit{v}}_1$ (m/day)	${\mathit{D}}_1$ (m${}^{\mathit{\beta }}$/day)	$\mathit{\beta }$	${\mathit{R}}_1$	${\mathit{v}}_2$ (m/day)	${\mathit{D}}_2$ (m${}^2$/day)	${\mathit{R}}_2$	${\mathit{R}}_3$	${\mathit{H}}_{\mathit{x}}\left(\mathit{n}\right)$
27									60
132	0.018	10	1.4	10	0.18	20	10	1	65
224									70
328									80

lists the best-fitting parameters of the PDDO-based and ADE models, respectively. The results show that the horizon $H_x$ of the PDDO-based model increases with time evolution, which is consistent with the experimental data.

The non-local pollutant transport that exists in various systems can be attributed to very different mechanisms. For example, in open water (e.g., rivers), the turbulence and the heterogeneous riverbed that will have an accelerated effect on the pollutant transport lead to non-local transport behavior. The early arrivals and trailing edges of pollutants (as shown in Figure 4) can be observed for highly permeable preferential paths in non-homogeneous media, especially in the media with fractures [18]. The ADE model does not consider the long-distance spatial correlation (i.e., the non-local nature in space), resulting in its inability to effectively capture the transport behavior of pollutants in non-homogeneous media. However, the F-ADE model assumes that pollutant particles can jump an infinite distance in one movement, which is inconsistent with the actual situation. The FADE model is a global domain model that may often overestimate the extent of the non-local transport of pollutants. More specifically, in non-homogeneous domains, highly developed spatial continuity causes the preferential flow path, producing high-velocity regions, and the pollutant tracers have a relatively high probability in the non-local jump. Due to the adjustability of the weight function and horizon zone, the PDDO-based model can describe pollutant non-local transport in various conditions. Considering the diversity of pollutant non-local transport mechanisms, the PDDO-based model is a generalized model that can describe the transport of pollutants in aquifers, the transport of pollutants in surface waters, etc.

4. Discussion

Diffusion Regime of the PDDO-Based Model

The theoretical issue of pollutant diffusion has long been central when using fluid mechanics to explore pollutant transport. A standard tool to explore the stochastic motion of Brownian tracers is the mean squared displacement (MSD) of the tracers. Similar to Brownian motion, this subsection aims to evaluate the PDDO-based model by comparing the MSD, which is the criteria for diffusion states in quantifying pollutant transport. The MSD is herein defined as follows [40]:

$$\begin{array}{c}\hfill R_{MSD}\left(t\right)=〈{\left(x\left(t\right)-〈x\left(t\right)〉\right)}^2〉\end{array}$$

(26)

Equation (26) gives an analytical form of the MSD of particles. Here, we use the following equation to calculate the MSD of particles in the framework of the PDDO-based model [31,40]:

$$\begin{array}{c}\hfill R_{MSD}\left(t\right)={\int }_{-\infty }^{+\infty }{x}^{2}c(x,t)dx\end{array}$$

(27)

where $c(x,t)$ is the solution of the model (9) in the case of $v=0$ with an instantaneous point source $c(x,0)=\delta \left(x\right)$.

Based on the analytical solution, (23) of model (9), and set $v=0$, yields

$$\begin{array}{cc}\hfill c(x,t)& =F^{-1}\left\{exp\left(\frac{1}{R}\left[D·{\widehat{g}}_N^2\left(k\right)-\beta \right]·t\right)\right\}\left(k\right)\hfill \\ \hfill \beta & =\underset{H_x}{\int }\left(D·g_N^2\left(\xi \right)\right)d\xi \hfill \end{array}$$

(28)

Considering a general form of ${\widehat{g}}_N^2\left(k\right)$ as follows

$$\begin{array}{c}\hfill {\widehat{g}}_N^2\left(k\right)=-C·{\left|k\right|}^{\alpha }\end{array}$$

(29)

Inserting Equation (29) into Equation (28), one obtains

$$\begin{array}{cc}\hfill c(x,t)& =F^{-1}\left\{exp\left(-C_1{\left|k\right|}^{\alpha }·t\right)\right\}\left(k\right)\hfill \end{array}$$

(30)

where $C_1=D·C/R$

First, we consider that $\alpha =2$, taking inverse Fourier transform of Equation (30) yields

$$\begin{array}{c}\hfill c(x,t)=\frac{1}{\sqrt{2C_{1}t}}exp\left(-\frac{x^2}{4C_{1}t}\right)\end{array}$$

(31)

Inserting Equation (31) into Equation (27), one obtains

$$\begin{array}{c}\hfill R_{MSD}\left(t\right)=2\sqrt{2\pi }{C}_{1}t\end{array}$$

(32)

Equation (32) gives that the MSD $R_{MSD}\left(t\right)$ is linear with time, and the Fickian diffusion regime is obtained in the case of $\alpha =2$.

Furthermore, considering that $1<\alpha <2$, Equation (30) is the characteristic function of a centered and symmetric Lévy distribution. Therefore, the $c(x,t)$ has the following power-law asymptotic:

$$\begin{array}{c}\hfill c(x,t)\sim \frac{C_{1}t}{{\left|x\right|}^{1+\alpha }}\end{array}$$

(33)

Due to Equation (33) being a Lévy distribution, the mean squared displacement diverges:

$$\begin{array}{c}\hfill 〈x^2\left(t\right)〉\to \infty \end{array}$$

(34)

Hence, the analytical result shows that the MSD is nonlinear with time in the case of $1<\alpha <2$, and the anomalous diffusion regime occurs. Based on the above analysis of the MSD of the PDDO-based model, the proposed non-local model can capture both normal (Fickian) and anomalous diffusion regimes of pollutants. The PDDO-based model bridges the local and regional scales by considering various PD functions.

It is also noteworthy that this study has three main limitations of the PDDO-based model (9) in quantifying pollutant non-local transport behavior. First, the weight function $w\left(\right|\xi \left|\right)$ is chosen as empirical, which should be further explored for various conditions. Second, the horizon $H_x$ is determined by fitting the experimental data, and the quantitative relationship between $H_x$ and media needs to be established in the future. Third, the non-local property leads to the calculation time of the PDDO-based model being much higher than that of the local model. Fast algorithms for solving the PDDO-based model need to be further investigated.

5. Conclusions

This study proposes a non-local model for quantifying pollutant transport in a heterogeneous environment. The non-local property of pollutant transport is encapsulated in the PDDO. Three conclusions can be drawn from this study.

(1) PDDO gives an efficient tool to describe the non-local transport behavior of pollutants, and the non-local transport feature of pollutants can be adjusted by considering different weight functions $w\left(\right|\xi \left|\right)$ and horizon $H_x$.

(2) The analyses of the MSD for the PDDO-based model give that the PDDO model can characterize both the normal (Fickian) and anomalous transport behavior of pollutants. The PDDO-based model bridges the local scale and region scale by considering various PD functions.

(3) Applications show that the PDDO-based model can efficiently characterize the observed non-local transport behavior of pollutants, and the non-local feature of pollutant transport enhances with time evolution.