Mathematical modeling serves as an indispensable tool for comprehending the dynamic mechanisms underlying the spread of infectious diseases and devising effective strategies to mitigate economic losses. A pivotal achievement in mathematical epidemiology lies in the formulation, prediction, and identification of disease control measures through rigorous mathematical expressions. Traditionally, epidemiological models have compartmentalized a constant population into three categories: susceptible, infectious, and recovered. The evolution of these models has seen the inclusion of various compartment classes, adapting to the variation of available descriptive information. In the budding stages of research, the discourse predominantly revolved around integer-order differential equations. However, recent strides in the field have witnessed the emergence of fractional-order mathematical modeling, playing a substantial role in addressing a diverse array of challenges in engineering, epidemiology, ecology, and biology. Fractional calculus, as a generalization where the conventional integer order gives way to a fractional order, provides a more varied approach . Notably, as the order approaches one, the solution of the fractional order converges to that of the integer order system. The fractional order system introduces two critical properties—the memory property and hereditary property—which are not adequately captured by any integer order system. These properties contribute to a more accurate representation of real-world phenomena. The fractional order system aligns well with theoretical descriptions of various mechanisms, offering a sophisticated understanding of dynamics . Over the last few decades, various types of fractional calculus, including Riemann–Liouville, Caputo, Atangana–Baleanu, Caputo–Fabrizio, Katugampola, Hadamard, etc., have been introduced extensively to scrutinize the dynamics of epidemic models. This ongoing exploration holds the promise of refining our capacity to model and comprehend intricate systems across diverse fields.

Subtrapaul, et., proposed the dynamical behavior of the COVID-19 epidemic SIQR model in terms of fractional order. S. S. Askar, and Dipankar Ghosh, studied the h-stability of the Caputo fractional order differential equation by using some fractional comparison principle and the Lyapunov direct method. Mukherjee D, and Mondal R, deal with the fractional order prey-predator model of the dynamical system for the reserved area and present the numerical analysis of hopf-bifurcation. Zizhen Zhang et al., formulated and analyzed the dynamics of the COVID-19 epidemic models with an isolated class of fractional order. Many researchers proposed many models in eco-epidemiology and analyzed the various stability properties, existence, and non-negativity of the solution . Numerical approximation of the proposed model has been done by the Rung-Kutta method, Taylor series method, Adams Bashforth and Adams Moulton method, Adomain and Laplace Adomian decomposition method, Homotopy method and generalization of homotopy method . In this work, the Laplace Adomain decomposition method is used to study the approximate solution of the proposed model. The Laplace Adomian decomposition method was introduced by Adomian in 1980. This is a very effective method to find the explicit and numerical solution to the physical problem of ordinary and partial differential equations with nonlinear initial and boundary conditions. It can also be used to study many deterministic, stochastic diffusion models with and without delay. There is no perturbation and linearization is required and also no extra memory is required to solve this problem which is the main advantage of the method .

Sensitivity analysis (SA) is increasingly popular in developing more efficient models, requiring reliable statistical and mathematical methods to enhance modeling accuracy. SA serves various purposes, including developing decisions or recommendations, communicating, understanding, or quantifying systems, and refining models. During model development, SA verifies validity or accuracy, simplifies and calibrates weak processes or missing data, and identifies important parameters for further studies. These applications highlight SA's versatility and importance in refining models and optimizing decision-making processes.

The main aim of this paper is to study the mathematical model of Caputo fractional order soil pollution. Soil is the main requirement for living and non-living organisms. In the case of contamination of soil as such oil spill, excess use of heavy metals, and inorganic fertilizer, the soils are polluted. All polluted soil contains compounds that include metals, inorganic iron salts (phosphates, carbonates, sulfates, nitrates), and organic compounds (lipids, DNA, fatty acids, alcohol). These compounds are mainly entered through soil microbial activity and decomposition of organisms (animals and plants) . The main soil pollution has an unfavorable two-way impact on food security, due to toxic contamination, which can reduce crop yields and be unsafe for consumption by humans and animals. Currently, the degradation of soil and land is affecting 3.2 million people which is 40% of the world’s population. It is necessary to treat the damaged soil and to maintain better soil management practices to limit soil pollution.

In this study, the work is carried out as follows: In Section 2, the mathematical model of the Caputo fractal fractional order SPTR model is presented. Section 3, analyses the existence, uniqueness, boundedness, and non-negative of the solution. In Section 4, we discuss the equilibrium and stability analysis. An adaptive control and sensitivity analysis are performed in Section 5. In Section 6, the numerical approximation of the solution of the model is done. In Section 7, the numerical simulation and discussion are given. The conclusion is given in Section 8.

The dynamic model for soil pollution comprises four distinct compartments: susceptible soil areas S(t), polluted soil areas (P), areas subject to necessary remedial treatments (T), and the recovered subclasses (R). Soil contamination typically arises from the presence of excessive heavy metals and plastics. Consequently, the application of both organic and inorganic remedies to the soil treatment (T) class becomes crucial, where certain areas may experience recovery while others might reintroduce contaminants to the susceptible soil due to treatment failures. In the contemporary landscape, there is a notable surge in the utilization of heavy metals in daily activities, leading to a recurrent cycle where previously recovered soil reverts to a susceptible state. In light of the above, the transmission dynamics of soil pollution under treatment measures are given in the compartmental representation illustrated in Figure 1. This model encapsulates the intricate interplay between susceptible and treated soil areas, offering insights into the complex dynamics of soil pollution and the efficacy of remedial measures.

The mathematical model of fractional order is presented as follows,

with the initial conditions,

where the parameter and their descriptions are as follows, $N$ is the total soil area with $N=S+P+T+R$

B – recruitment rate of the soil,

$\alpha$ - the pollution rate of the soil,

$\beta$- recovered re-enter to the susceptible rate due to failure of the treatment,

$\gamma$ – natural recovered rate,

${\alpha }_1$ - remediation rate,

${\gamma }_1$- recovered rate due to treatment,

$\delta$- recovered re-enter to the susceptible class,

$\mu$- death rate of all the classes.

For the fractional order system, the existence and uniqueness of solutions in the region $R_M\times \left(0,t\right)$, where $R_M=\left\{\left(S,P,T,R\right)\in R_{+}^4:\max \left(\left|S\right|,\left|P\right|,\left|T\right|,\left|R\right|\le \eta \right)\right\}$ is considered, where $\eta$ is some positive constant.

Theorem 3.1 If${\mathit{ Y}}_0=\left(S_0, P_0,T_0,R_0\right)\in R_M,$then for each $Y_0$ there exists a unique solution $Y\left(t\right)\in R_M$ of the system (1) with initial condition $Y_0$ for all $t\ge 0$.

Proof: Let us consider the approach of and take $R\left(Y\right)=(R_1\left(Y\right),R_2\left(Y\right),R_3\left(Y\right),R_4(Y\left)\right)$,

The system of equations (3) is reduced to,

where, $A=\gamma {+\alpha }_1+\mu$, $\mathit{ B}=\beta +{\gamma }_1+\mu$, $C=\delta +\mu$. For any$,\check{Y}\in R_M$,

where $H=\max \{ \left(2\mathit{\alpha \eta }+\mu \right),\left(2\left(\mathit{\alpha \eta }+{\alpha }_1+\gamma \right)+\mu \right),\mathit{ }2\left(\beta +{\gamma }_1\right)+\mu ,\mathit{ }2\delta +\mu \}$

Hence $R\left(Y\right)$ satisfies the Lipshitz condition and so the existence and uniqueness of the solution of fractional order system (1) with the initial conditions are confirmed.

Theorem 3.2 The system of the fractional order of (1) in $R_M$ is non-negative and uniformly bounded.

Proof: The approach used by is followed. Define the function $W\left(t\right)=S\left(t\right)+P\left(t\right)+T\left(t\right)+R\left(t\right)$ and

Applying the standard comparison theorem for the fractional order given in , we have,

where $E_{\zeta }$ is the Mittag Leffler function. By taking $t\to \infty$, $0\le W\left(t\right)\le \frac{B}{\mu }$. Therefore, the solution is uniformly bounded in the region $R=\left\{\left(S,P,T,R\right)\in R_{+}^4,W\left(t\right)\le \frac{B}{\mu }+ϵ,\mathit{ ϵ}>0 \right\}.$

Now the non - non-negativity of the solution is studied for the fractional order system (1). In (1)

Hence the solution of the system is non-negative.

There are two equilibrium points exists for the model that are found by equating time derivatives in system (1) to zero as follows,

The two steady states of the model (1) are given by,

Pollution free equilibrium point: $P_0\left(\frac{B}{\mu },0,0,0\right)$

Pollution extinct equilibrium point: $P_{*}(S^{*},P^{*},T^{*},R^{*})$

where $S^{*}=\frac{\gamma +{\alpha }_1+\mu }{\alpha }$, $T^{*}=\frac{{\alpha }_1{P}^{*}}{\beta +{\gamma }_1+\mu }$, $R^{* }=\frac{1}{\delta +\mu }\left[\gamma +\frac{{\gamma }_1{\alpha }_1}{\beta +{\gamma }_1+\mu }\right]P^{*}$

$P_{*}$ is positive, when (i) $R_0<1$ and $\gamma +{\alpha }_1+\mu <\left(\frac{\beta {\alpha }_1}{\beta +{\gamma }_1+\mu }+\frac{1}{\delta }\left(\frac{1}{\delta +\mu }\left(\gamma +\frac{{\gamma }_1{\alpha }_1}{\beta +{\gamma }_1+\mu }\right)\right)\right)$ (or)

(ii) $R_0>1$ and $\gamma +{\alpha }_1+\mu >\left(\frac{\beta {\alpha }_1}{\beta +{\gamma }_1+\mu }+\frac{1}{\delta }\left(\frac{1}{\delta +\mu }\left(\gamma +\frac{{\gamma }_1{\alpha }_1}{\beta +{\gamma }_1+\mu }\right)\right)\right)$

Basic Reproduction Number (BRN):

The next generation matrix (NGM) is used to find the BRN of the system (1)

Consider,${}_t{}^{\mathit{CFF}}{\mathbb{D}}^{\mathit{\zeta },\wp }A\left(t\right)=F\left(Y\right)-V\left(Y\right)$, where$\mathit{ F}\left(Y\right)=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=\left[\begin{array}{c}\mathit{\alpha SP}\\ 0\\ 0\end{array}\right]$

where $A=-(\gamma +{\alpha }_1+\mu ),\mathit{ D}=-(\beta +{\gamma }_1+\mu ) ,\mathit{ C}=-(\delta +\mu ).$It follows that spectral radius of the matrix is $\rho \left(F{V}^{-1}\right)=\max \left({\lambda }_{1,2,3}\right)$ at pollution free equilibrium point $P_0$. Therefore, the basic reproduction number is,

Theorem: 4.1 If all the eigen values have negative real part then the pollution free equilibrium point $P_0(S_0,0,0,0)$ of the system (1) is locally asymptotically stable.

Proof: From the Matignon condition, we observe that the pollution free equilibrium point is locally asymptotically stable if and only if all the eigenvalues ${\lambda }_i$ of the Jacobian $\left(J\right(P_0\left)\right)$ satisfy $|\arg ({\lambda }_i\left)\right|>\frac{\mathit{a\pi }}{2}.$ The Jacobian matrix of the system (1) is

Here are the eigenvalues, ${\lambda }_1=-\mu ,\mathit{ }{\lambda }_2=-\left(\delta +\mu \right), {\mathit{ \lambda }}_{3,4}=\left[\begin{array}{cc}-\left(\beta +\gamma +\mu \right)& 0\\ {\gamma }_1& -(\delta +\mu )\end{array}\right]$. Therefore, the eigenvalues ${\lambda }_3=-\left(\beta +\gamma +\mu \right), {\mathit{ \lambda }}_4=-(\delta +\mu )$. It can be observed that all the eigenvalues of the Jacobian matrix are negative. Hence by using the Matignon condition, we observe that the pollution-free equilibrium is locally asymptotically stable.

Theorem 4.2: The pollution extinction equilibrium point${\mathit{ P}}_{*}\left(S^{*},P^{*},T^{*},R^{*}\right)$ is globally asymptotically stable if and only if $R_0>1$ and $\frac{\mu }{S}+\frac{{\beta }^{2}T{T}^{*}}{4\alpha }+\frac{{\delta }^{2}R{R}^{*}}{4(\gamma +{\gamma }_1)}\le \frac{P^{*}}{S}.$

Proof: Consider the following nonlinear Lyapunov function,

Take the Lyapunov fractional derivative on both sides,

After some simplification we obtain,

Finally, we conclude that, ${}_t{}^{\mathit{CFF}}{\mathbb{D}}^{\mathit{\zeta },\wp }W\left(t\right)\le 0$, if and only if $\frac{\mu }{S}+\frac{{\beta }^{2}T{T}^{*}}{4\alpha }+\frac{{\delta }^{2}R{R}^{*}}{4(\gamma +{\gamma }_1)}\le \frac{P^{*}}{S}.$

The first derivative of a Lyapunov function is instrumental in evaluating the global stability of endemic equilibrium locations. It provides essential insights that can be further enriched by second derivative analysis. While the second derivative informs us about curvature through its sign, the first derivative offers valuable information on disease propagation. Examining the first derivative gives us insights into how the system evolves and spreads over time. Complemented by the second derivative examination, this analysis contributes to a comprehensive understanding of the dynamics and stability of endemic equilibria in disease models.

Taking the second derivative of (9),

where,

Now, substituting the values of the derivative $\dot{S},\dot{P},\dot{T},\dot{R}$ and rearranging the terms of positive and negative groups, so that we have,

Therefore, we see that,

i. If ${\mathcal{H}}_1>{\mathcal{H}}_2$ then $\frac{d^{2}W\left(t\right)}{\mathit{dt}}>0,$

ii. If ${\mathcal{H}}_1<{\mathcal{H}}_2$ then $\frac{d^{2}W\left(t\right)}{\mathit{dt}}<0,$

iii. If ${\mathcal{H}}_1={\mathcal{H}}_2$ then $\frac{d^{2}W\left(t\right)}{\mathit{dt}}=0,$

This interpretation explains the model's global behavior with respect to Lyapunov second order derivative.

Lemma 4.1: The given model (1) has no periodic orbits.

Proof: The proof has carried with a Dulac plus Poincare Bendixson theorem as follows:

$h\left(S,P\right)=\frac{1}{\mathit{SP}}$ , where all the parameters are positive.

Hence, from the Dulac principle, we conclude that there is no periodic solution for the given region R.

The numerical approximate solution of a system (1) is obtained using Laplace adomain decomposition method,

Applying Laplace transform on both sides of equation (1),

Substitute the initial condition (2) in the above equations,

Now,

Consider the infinite series,

Also consider the decomposition of nonlinear terms as follows,

where the decomposition polynomial $A_k$ is defined as follows,

The first three polynomials are,

Substituting above equations in equation (11),

Comparing both sides of equation of (14),