Assessing contamination sources and environmental hazards for potentially toxic elements and organic compounds in the soils of a heavily anthropized area: the case study of the Acerra plain (Southern Italy)

1.   Introduction

In a letter to L'Hospital in 1695, the famous mathematician Leibniz introduced the concept of fractional derivative for the first time. Fractional calculus is concerned with non-integral differential and integral operators. The integer-order differential operator is a local operator, but the fractional-order differential operator is non-local, indicating that a system's next state is determined not only by its current state but also by all of its past states. It is more realistic, which is one of the main reasons for the popularity of fractional calculus. Fractional calculus was found to be more suitable for modelling real-world problems than classical calculus. The theory of fractional calculus provides an effective and systematic interpretation of nature's reality. The following are some basic works on fractional calculus on various topics, such as Podlubny ^[1], Caputo ^[2], Kiryakova ^[3], Jafari and Seifi ^[4,5], Momani and Shawagfeh ^[6], Oldham and Spanier ^[7], Diethelm et al. ^[8], Miller and Ross ^[9], Kemple and Beyer ^[10], Kilbas and Trujillo ^[11].

Fractional differential equations (FDEs) are widely widely used in various fields of science. FDEs have attracted a lot of attention in the last few years because of their diverse uses in physics and engineering. Due to its proven usefulness in a wide range of very diverse disciplines of science and engineering, fractional partial differential equations (PDEs) have gained importance and reputation among FDEs in recent years. For example, fractional derivatives in a fluid-dynamic traffic model can be utilised to overcome the insufficiency produced by the assumption of a continuous flow of traffic. Nonlinear FPDE solutions are of tremendous interest in both mathematics and practical applications ^{[12,13,14,15,16,17,18,19,20,21]}. The world's most important processes are represented by nonlinear equations. Manipulation of nonlinear processes is essential in physics, applied mathematics, and engineering problems. The significance of finding the exact solution to nonlinear partial differential equations is still a prominent issue in physics and applied mathematics, requiring implemented using different techniques to determine new approximate or exact solutions ^{[22,23,24,25]}.

In order to obtain explicit solutions for nonlinear equations of integer order, various techniques have been used. To solve FPDEs, however, only a less numerous approaches are used. such as Laplace variational iteration method (LVIM) ^[26], iterative Laplace transform method (ILTM) ^[27], optimal Homotopy asymptotic method (OHAM)^[28], approximate-analytical method (AAM) ^[29], reduced differential transform method (RDTM) ^[30], Laplace-Adomian decomposition method (LADM) ^[31], natural transform decomposition method (NTDM) ^[32], Elzaki transform decomposition method (ETDM) ^[33,34,35], new meshfree technique (NMT) ^{[36,37,38,39,40]}, Homotopy analysis method (HAM) ^[41], generalized exponential rational function method (GERFM) ^[42,43] and much more ^{[44,45,46,47,48,49,50,51]}. The goal of this paper is to implement the natural decomposition approach to solve modified Boussinesq equations and approximate long wave equations having fractional-order. Natural decomposition methods avoid roundoff errors by not requiring prescriptive assumptions, linearization, discretization, or perturbation.

Whitham ^[52], Broer ^[53], and Kaup ^[54] discovered the Whitham-Broer-Kaup (WBK) equations in the twentieth century, which describe the propagation of shallow water waves with various dispersion relations. Consider the fractional order coupled WBK equations ^[55]:

where $v = v(\xi, \tau)$ is the height that deviates from the liquid's equilibrium position and $u = u(\xi, \tau)$ is the horizontal velocity. The order of the time-fractional derivative is in this case. Furthermore, a and b are constants that represent different diffusion strengths, therefore Eq (1.1) becomes a modified Boussinesq equation if $a = 1$ and $b = 0$. Similarly, the system denotes the standard long wave equation for $a = 0$ and $b = 1.$ These equations are used in ocean and coastal engineering to explain the propagation of waves in dissipative and nonlinear media. They are recommended for situations involving water leakage in porous subsurface stratum and are frequently used in hydrodynamics. Furthermore, Eq (1.1) serves as the foundation for a number of models that represent the unconfined subsurface, such as drainage and groundwater problems ^[56,57].

The following is a summary of the paper's structure. The essential definitions related to fractional calculus are briefly discussed in Section 2. The NTDM is a method implemented to solve coupled WBK equations having fractional order with non-singular definitions are given in Section 3. In Section 4, we discussed the uniqueness and convergence results. Two examples of fractional-order modified Boussinesq and approximate long wave equations are given in Section 5 to validate the approaches. The article's brief conclusions are offered in Section 6.

2.   Basic preliminaries

In the literature, there are various fractional derivative definitions; for additional details, see ^[58,59,60]. This section includes Caputo, Riemann-Liouville, Atangana-Baleanu and Caputo-Fabrizio definitions for the benefit of the readers.

Definition 2.1. For a function $h\in C_{v}, v\geq-1$, the Riemann-Liouville fractional integral operator is given as ^[61]

Definition 2.2. The fractional derivative $h(\eta)$ in Caputo manner is given as ^[61]

for ${m}-1 < {{\wp}}\leq {m}, \ \ {m}\in N, \ \ \eta > 0, h\in C^{{m}}_{v}, v\geq-1$.

Definition 2.3. Caputo-Fabrizio fractional derivative for $h(\eta)$ is given as ^[61]

where $0 < {{\wp}} < 1$ and $F({{\wp}})$ is a normalization function with $F(0) = F(1) = 1$.

Definition 2.4. The Atangana-Baleanu Caputo derivative of fractional-order for $h(\eta)$ is defined as ^[61]

where $0 < {{\wp}} < 1$, where $B(\wp)$ is Normalization function and $E_{{{\wp}}}(z) = \sum_{m = 0}^{\infty}\frac{z^{m}}{\Gamma({{\wp}}^{m}+1)}$ is the Mittag-Leffler function.

Definition 2.5. For a function $u(\Im)$, the Natural transformation is defined as

Natural transformation of $u(\Im)$ for $\Im\in(0, \infty)$ is given as

here $H(\Im)$ is the Heaviside function.

Definition 2.6. For a function $\mathcal{U}({s}, {{\vartheta}})$, the Natural inverse transformation is given as

Lemma 2.1. If linearity property having Natural transformation for $u_{1}(\Im)$ is $u_{1}({s}, {{\vartheta}})$ and $u_{2}(\Im)$ is $u_{2}({s}, {{\vartheta}})$, then

where $c_{1}$ and $c_{2}$ are constants.

Lemma 2.2. If Natural inverse transformation of $\mathcal{U}_{1}({s}, {{\vartheta}})$ and $\mathcal{U}_{2}({s}, {{\vartheta}})$ are $u_{1}(\Im)$ and $u_{2}(\Im)$ respectivelythen

where $c_{1}$ and $c_{2}$ are constants.

Definition 2.7. The Natural transformation of Caputo operator $D_{\Im}^{{{\wp}}}u(\Im)$ is define as ^[61]

Definition 2.8. In terms of Caputo-Fabrizio, the Natural transformation of $D_{\Im}^{{{\wp}}}u(\Im)$ is given as ^[61]

Definition 2.9. In terms of Atangana-Baleanu Caputo derivative, the Natural transformation of $D_{\Im}^{{{\wp}}}u(\Im)$ is defined as ^[61]

3.   Methodology

This section introduces a general numerical methodology based on the Natural transform for the following equation.

with initial source

here $\mathcal{L}$ is linear, $\mathbb{N}$ is nonlinear and $h(\xi, \Im)$ represents source term.

3.1.   Case Ⅰ $(NTDM_{CF}):$

We get this by applying Caputo-Fabrizio fractional derivative with the aid of Natural transformation of Eq (3.1),

where

We may express Eq (3.3) as, using the inverse Natural transformation

$\mathbb{N}(\nu(\xi, \Im))$ can be decomposed into

where the Adomian polynomials is represented by $A_{\Im}$.We assume that the numerical solution of Eq (3.1) exists

By putting Eqs (3.6) and (3.7) into (3.5), we get

From (3.8), we obtain

we obtain the $NTDM_{CF}$ solution of (3.1) by putting (3.9) into (3.7),

3.2.   Case Ⅱ $(NTDM_{ABC}):$

We get this by applying Atangana-Baleanu derivative with the aid of Natural transformation of Eq (3.1),

where

On applying inverse Natural transformation (2.7), we express (3.11) as,

$\mathbb{N}(\nu(\xi, \Im))$ can be decomposed into

where the Adomian polynomials^[62,63] is represented by $A_{\Im}$.We assume that the numerical solution of Eq (3.1) exists

By putting Eqs (3.14) and (3.15) into (3.13), we get

From (3.8), we get

we obtain the $NTDM_{ABC}$ solution of (3.1) by putting (3.17) into (3.15),

4.   Convergence analysis

Here we discuss uniqueness and convergence and of the $NTDM_{CF}$ and $NTDM_{ABC}$.

Theorem 4.1. The result of (3.1) is unique for $NTDM_{CF}$ when $0 < (\lambda_{1}+\lambda_{2})(1-{{\wp}}+{{\wp}} \Im) < 1.$

Proof. Let $H = (C[J], ||.||)$ with the norm $||\phi(\Im)|| = {max}_{\Im\in J}|\phi(\Im)|$ is Banach space, $\forall$ continuous function on $J.$ Let $I:H\rightarrow H$ is a non-linear mapping, where

Suppose that $|\mathcal{L}(\nu)-\mathcal{L}(\nu^{*})| < \lambda_{1}|\nu-\nu^{*}|$ and $|\mathbb{N}(\nu)-\mathbb{N}(\nu^{*})| < \lambda_{2}|\nu-\nu^{*}|$, where $\nu: = \nu(\zeta, \Im)$ and $\nu^{*}: = \nu^{*}(\zeta, \Im)$ are are two different function values and $\lambda_{1}$, $\lambda_{2}$ are Lipschitz constants.

I is contraction as $0 < (\lambda_{1}+\lambda_{2})(1-{{\wp}}+{{\wp}} \Im) < 1$. From Banach fixed point theorem the result of (3.1) is unique.

Theorem 4.2. The result of (3.1) is unique for $NTDM_{ABC}$ when $0 < (\lambda_{1}+\lambda_{2})(1-{{\wp}}+{{\wp}}\frac{\Im^{mu}}{\Gamma(mu+1)}) < 1.$

Proof. Let $H = (C[J], ||.||)$ with the norm $||\phi(\Im)|| = {max}_{\Im\in J}|\phi(\Im)|$be the Banach space, $\forall$ continuous function on $J.$ Let $I:H\rightarrow H$ is a non-linear mapping, where

Suppose that $|\mathcal{L}(\nu)-\mathcal{L}(\nu^{*})| < \lambda_{1}|\nu-\nu^{*}|$ and $|\mathbb{N}(\nu)-\mathbb{N}(\nu^{*})| < \lambda_{2}|\nu-\nu^{*}|$, where $\nu: = \nu(\zeta, \Im)$ and $\nu^{*}: = \nu^{*}(\zeta, \Im)$ are are two different function values and $\lambda_{1}$, $\lambda_{2}$ are Lipschitz constants.

$I$ is contraction as $0 < (\lambda_{1}+\lambda_{2})(1-{{\wp}}+{{\wp}}\frac{\Im^{{{\wp}}}}{\Gamma{{{\wp}}+1}}) < 1$. From Banach fixed point theorem the result of (3.1) is unique.

Theorem 4.3. The $NTDM_{CF}$ result of (3.1) is convergent.

Proof. Let $\nu_{m} = \sum_{r = 0}^{m}\nu_{r}(\zeta, \Im)$. To show that $\nu_{m}$ is a Cauchy sequence in H. Let,

Let $m = n+1$, then

where $\lambda = (\lambda_{1}+\lambda_{2})(1-{{\wp}}+{{\wp}} \Im)$. Similarly, we have

As $0 < \lambda < 1$, we get $1-\lambda^{m-n} < 1$. Therefore,

Since $||\nu_{1}|| < \infty, \ \ ||\nu_{m}-\nu_{n}||\rightarrow 0$ when $n \rightarrow \infty$. As a result, $\nu_m$ is a Cauchy sequence in H, implying that the series $\nu_m$ is convergent.

Theorem 4.4. The $NTDM_{ABC}$ result of (3.1) is convergent.

Proof. Let $\nu_{m} = \sum_{r = 0}^{m}\nu_{r}(\zeta, \Im)$. To show that $\nu_{m}$ is a Cauchy sequence in H. Let,

Let $m = n+1$, then

where $\lambda = (\lambda_{1}+\lambda_{2})(1-{{\wp}}+{{\wp}}\frac{\Im^{{{\wp}}}}{\Gamma({{\wp}}+1)})$. Similarly, we have

As $0 < \lambda < 1$, we get $1-\lambda^{m-n} < 1$. Therefore,

Since $||\nu_{1}|| < \infty, \ \ ||\nu_{m}-\nu_{n}||\rightarrow 0$ when $n \rightarrow \infty$. As a result, $\nu_m$ is a Cauchy sequence in H, implying that the series $\nu_m$ is convergent.

5.   Numerical examples

In this section we investigate the analytical solution for a few problems of fractional order coupled WBK equations. We chose these equations because they contain closed form solutions and are well-known methods for analysing results in the literature.

5.1.   Example 1

Consider the modified Boussinesq (MB) equations of fractional order

with initial source

By applying the Natural transform to Eq (5.1), we have

Define the non-linear operator as

The above equation is reduced to by simplifying it

We get by applying inverse $NT$ to Eq (5.5)

5.1.1.   Now we implementing $NTDM_{CF}$

Assume that the unknown functions $u(\xi, \Im)$ and $v(\xi, \Im)$ yield the infinite series result as follows

The Adomian polynomials represent the non-linear terms and are denoted by $uu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal A_{{l}}$, $uv_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal B_{{l}}$ and $vu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal C_{{l}}$. We can write the Eq (5.6) as a result of applying these terms,

We can write as follows by comparing both sides of Eq (5.8)

The remaining components of the Natural decomposition technique result $u_{{l}}$ and $v_{{l}}for ({l}\geq3)$ can be performed easily. The series form result is calculated as follows:

5.1.2.   Now we implementing $NTDM_{ABC}$

Assume that the unknown functions $u({\xi}, {\Im})$ and $v({\xi}, {\Im})$ yield the infinite series result as follows

The Adomian polynomials represent the non-linear terms and are denoted by $uu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal A_{{l}}$, $uv_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal B_{{l}}$ and $vu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal C_{{l}}$. We can write the Eq (5.6) as a result of applying these terms,

We can write as follows by comparing both sides of Eq (5.14)

The remaining components of the natural decomposition technique result $u_{{l}}$ and $v_{{l}} for ({l}\geq3)$ can be performed easily. The series form result is calculated as follows:

For Eq (5.1), the exact result is obtained at $\wp = 1$,

5.2.   Example 2

Consider the approximate long wave (ALW) equations with arbitrary order

with initial source

By applying the Natural transform to Eq (5.20), we have

Define the non-linear operator as

The above equation is reduced to by simplifying it

We get by applying inverse $NT$ to Eq (5.24)

5.2.1.   Now we apply $NTDM_{CF}$

Assume that the unknown functions $u(\xi, \Im)$ and $v(\xi, \Im)$ yield the infinite series result as follows

The Adomian polynomials represent the non-linear terms and are denoted by $uu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal A_{{l}}$, $uv_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal B_{{l}}$ and $vu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal C_{{l}}$. We can write the Eq (5.25) as a result of applying these terms,

We can write as follows by comparing both sides of Eq (5.27)

The remaining components of the Natural decomposition technique result $u_{{l}}$ and $v_{{l}}for ({l}\geq3)$ can be performed easily. The series form result is calculated as follows:

5.2.2.   Now we apply $NTDM_{ABC}$

Assume that the unknown functions $u(\xi, \Im)$ and $v(\xi, \Im)$ yield the infinite series result as follows

The Adomian polynomials represent the non-linear terms and are denoted by $uu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal A_{{l}}$, $uv_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal B_{{l}}$ and $vu_{{\xi}} = \sum_{{l} = 0}^{\infty}\mathcal C_{{l}}$. We can write the Eq (5.25) as a result of applying these terms,

We can write as follows by comparing both sides of Eq (5.30)

The remaining components of the Natural decomposition technique result $u_{{l}}$ and $v_{{l}}({l}\geq3)$ can be performed easily. The series form result is calculated as follows:

For Eq (5.20), the exact result is obtained at $\wp = 1$

5.3.   Numerical results and discussion

We used two unique methods to investigate the numerical solution of systems of coupled modified Boussinesq and approximate long wave equations with fractional order in this work. Through Maple, you can find numerical data for the system of coupled modified Boussinesq and approximate long wave equations for any order at various values of space and time variables. For the system in Problem 1, we construct numerical simulations at various values of $\xi$ and $\Im$ in Tables 1 and 2. Tables 3 and 4 show a numerical comparison of the variational iteration method, Adomian decomposition method and Natural decomposition method in terms of absolute error for Eq (5.1).

The outcomes of a calculations for the coupled system considered in Problem 2 are shown in Tables 5 and 6. Similarly, in Tables 7 and 8, we compare the Adomian decomposition method, natural decomposition method and variational iteration method solutions to Eq (5.20).

On the basis of the information in the above tables, we may say that the Natural decomposition approach produces more accurate results.

Figure 1 illustrates the behaviour of the exact and Natural decomposition technique result for $u(\xi, \Im)$ of Problem 1, while Figure 2 displays the behavior of the analytical result at different fractional-orders of $\wp$. Figure 3 illustrates the solution for $u(\xi, \Im)$ in various fractional-orders. Figure 4 demonstrates the behaviour of the exact and analytical solutions for $v(\xi, \Im)$ and Figure 5 shows the absolute error for Example 1.

Figure 6 shows the behaviour of the exact and Natural decomposition technique results for $u(\xi, \Im)$ of Problem 2, whereas Figure 7 displays the nature of the analytical solution at various fractional-orders of $\wp$. Figure 8 shows the solution graph for Problem 2 at various fractional-orders of $u(\xi, \Im)$. Figure 9 shows the behavior of the exact and analytical solution for $v(\xi, \Im)$ whereas Figure 10 shows the absolute error of Problem 2. We draw the given graphs within domain $-100\leq\xi\leq100$ having $c = 10,$ $\omega = 0.005$ and $\ell = 0.15.$

6.   Conclusions

The natural decomposition technique is used to solve the coupled modified Boussinesq and approximate long wave equations in this paper. Two examples are solved to demonstrate and validate the effectiveness of the suggested technique. In comparison to existing analytical approaches for determining approximate solutions of nonlinear coupled fractional partial differential equations, the present method is efficient and straightforward. The derived results have been shown in graphical and tabular form. The proposed method generates a sequence of results in the form of a recurrence relation having greater accuracy and less calculations. In terms of absolute error, calculations were performed for both fractional coupled systems. At $\wp = 1$, a number of computational solutions are compared to well-known analytical techniques and the exact results. The benefits of the current methods include less calculations and greater accuracy. Moreover, the proposed method is shown to be simple and effective, and it may be applied to solve various fractional-order differential equation systems.

Acknowledgements

The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by Grant Code: (22UQU4310396DSR07).

Conflict of interest

The authors declare that they have no competing interests.