Numerical Schemes for Fractional Energy Balance Model of Climate Change with Diffusion Effects

This study aims to propose numerical schemes for fractional time discretization of partial differential equations (PDEs). The scheme is comprised of two stages. Using von Neumann stability analysis, we ensure the robustness of the scheme. The energy balance model for climate change is modified by adding source terms. The local stability analysis of the model is presented. Also, the fractional model in the form of PDEs with the effect of diffusion is given and solved by applying the proposed scheme. The proposed scheme is compared with the existing scheme, which shows a faster convergence of the presented scheme than the existing one. The effects of feedback, deep ocean heat uptake, and heat source parameters on global mean surface and deep ocean temperatures are displayed in graphs. The current study is cemented by the fact-based popular approximations of the surveys and modeling techniques, which have been the focus of several researchers for thousands of years.

about a 3 mm rise in sea level has been observed, and half of this rise is due to the expansion of water by heat [1].The rise in the earth's temperature is expected to cause an energy imbalance because it absorbs more solar radiation than it emits thermal radiation.In this regard, greenhouse gases play a significant role by forming a thick layer around the atmosphere and preventing thermal radiation's movement across it [7].Global warming increases the earth's temperature and increases the likelihood of volcanic eruptions, heat stroke, and severe radiation disproportion.Therefore, the consequences of global warming go beyond merely raising the temperature [8].
A luxuriant and mechanical lifestyle followed by clearing land for residential purposes has added more CO2 to the atmosphere.As a result, the whole earth acted as glass and raised the chances of becoming a greenhouse.This greenhouse gas formed an additional layer across the earth's surface and hindered the sun's rays from escaping it, increasing the temperature globally and leading to heatstroke or a heatwave [9][10][11].These temperature variations affect aquatic life and their movements, as evident from the shifted records of phytoplankton towards the North Sea.In contrast, this migration has changed the oceanic weather poleward [12,13].The motion of phytoplankton from the Atlantic to the cold water of green land has been variable so far; therefore, the ecologist and weather investigators collectively worked to control the migration and weather shift.In this regard, mathematicians have offered their services to calculate the mass change of phytoplankton through modeling and provide a suitable solution to prevent their extinction due to an unfavorable environment [14][15][16][17][18][19][20].
Later on, modern and fractional calculus have taken over the traditional one and can effectively elaborate on the present, past, and future [21][22][23][24][25].These studies revealed the practical conversion of real-world problems through the differential equation, which describes the whole density and spectrum lying between two integer orders.Modern and non-local fractional-order derivatives have ruled over conventional-order derivatives in fractional calculus.Firstly, Riemann and Liouville have proposed the first-ever definition of a fractional order derivative [26].Secondly, Caputo modified the conventional derivative by introducing a local singular kernel named the Caputo derivative.After that, in 2015, Caputo and Fabrizio further studied and manipulated this derivative by replacing the singular kernel, i.e., including the Mittag-Leffler, which is non-singular [27].Finally, in 2016, further modification of the Caputo and Fabrizio derivatives was done by building non-local and non-singular kernels by adding the Mittag-Leffler function [28].Nowadays, this fractional derivative has its worth in solving real-world problems [29,30].The Caputo derivative has the unique characteristic of having a single kernel.It grows at both ends of the transmission spectrum and exhibits a cross-over pattern for the amplitude of its mean square.This property enables the modeling of competition among living species for their survival.Climate change and pandemics can be modeled [31,32].A second-order numerical approximation for Caputo-Fabrizio fractional derivative has been given in Liu et al. [33] at node  + 1

2
. It was mentioned ( 2 ) operations and () Stages were required for the non-local property of the Caputo-Fabrizio, where the number of divided intervals was denoted by .Another fast algorithm based on a novel approximation technique was developed, which consumed () operations and (1) for memory stages.Another numerical approach, the local discontinuous Galerkin method [34], was employed with a new Caputo-Fabrizio fractional derivative for solving fractal mobile/immobile transport equations.The stability of the local discontinuous Galerkin method was proved, and prior error estimates were also derived.
The current study will provide great insight into preventing climatic change by using Caputo's fractional derivative, which is more efficient than the traditional integer model.The energy balance model is modified using the source term of the non-linear type.Also, the effect of the diffusion term describes the spread of temperature over the spatial variable.The partial differential equation version is more flexible than the standard model.The proposed scheme is capable of solving the modified model.The scheme is finite difference discretization for the temporal term in the parabolic equation.Results are supported through numerical simulation.
Predicting the climate system's evolution in response to natural and human-caused perturbations is impossible without resorting to mathematical models and numerical methods.Different types of mathematical models and numerical methods have played an important role for several decades.In Soldatenko et al. [35], for studying the earth's climate change, some methods and mathematical models were presented with the use of sensitivity analysis.The variation of climate was studied using stochastic models, and essential features of this modeling have been discussed.The application of a low-order energy balance model was explored for studying variations in climate.The influences of variations in feedback and the climate system's inertia were estimated using the sensitivity analysis approach.A simple radiative balance climate model has been presented in North [36] that includes constant homogeneous cloudiness, zonal averaging, an ice feedback mechanism, and ordinary diffusive thermal heat transfer.The simplest model with one parameter was solved analytically, and the solution was found in hypergeometric functions.The ice sheet latitude in terms of the solar constant was studied using this solution.It was proved that the presented climate and ice-covered earth were stable, which was done by stability analysis of the equilibrium solutions.
Mathematical modeling of the earth's climate is a powerful tool for studying climate variations in the future.Mathematical models construct sophisticated models in which a dynamical system represents the simulated processes and ordinary or partial differential equations are used to describe them.There are four climate models: a simple climate model, a general circulation model of the atmosphere based on sea ice and the upper mixed ocean layer, an intermediate complexity model, and a fully coupled atmosphere-ocean general circulation model.For accurate and reliable results, complex models can be proposed.Based on some factors, a suitable climate model can be chosen.Among the mathematical models of climate change, the energy balance model (EBM) is quite effective and simple.Budyko [37] and Sellers [38] introduced these models about 50 years ago.
The model equations for the energy balance model are expressed as: where  denotes the global mean surface temperature and   shows global mean ocean temperature,  represents the heat capacity of the upper box, which is consisted of oceanic and atmosphere surface mixed layer;  represents the climate feedback parameter which is a proportionality coefficient between the global mean surface temperature (GMST) and radiative response,   represents the effective heat capacity of the deep ocean,  is used for the deep-ocean heat uptake parameter, and  is the forcing term.
Since Equations 1 and 2 are ordinary differential equations that only use the information of time variables, classical ordinary time derivatives are also used in these equations, which are replaced by a more general fractional derivative with the diffusive spatial term.The diffusive term is used to describe the diffusion temperature over a spatial variable.The proposed fractional model can be expressed as: (3) where Equations 3 and 4 are fractional partial differential equations having diffusive terms and a source term in Equation 3. At this stage, a local stability analysis for the corresponding ordinary differential equations is presented for Equations 3 and 4 using  1 = 0 =  2 .
To find the equilibrium points, set 0 = − − ( −   ) +  2  (5) Solving Equations 5 and 6 give the following points of equilibrium: where  and   denotes, respectively first and second coordinates in  0 and  1 .
Since Equations 5 and 6 are non-linear, their Jacobin, if found, is expressed as: At the equilibrium point  1 The Jacobin can be expressed as: The characteristic polynomial () can be expressed as: where  denotes the eigenvalues of |  1 .For the system to be a stable eigenvalue of |  1 should be non-negative.To impose the condition of getting negative eigen values, the Routh Hurwitz criterion will be applied.According to this criterion for second-degree polynomials, the coefficients of  and the constant term of () must be negative.Therefore, 2 −  < 0 and − < 0 Since the constant term in Equation 9 is already negative, so the condition of getting a negative eigenvalue is  < 2.So, when this mentioned condition is satisfied, the system of Equations 3 and 4 using  1 = 0 =  2 will be stable.Equations 3 and 4 will be solved numerically using a numerical fractional scheme.

Basic Concepts:
Before starting the working procedure of numerical schemes for fractional differential equations, some basic literature related to the fractional derivative is given.

Definition:
The Caputo fractional derivative is defined as [39]: Definition: For a function , the fractional integral operator is defined as: Lemma: If  − 1 <  ≤ , then the following properties hold: Since the numerical scheme for the fractional differential equation is already given in Soldatenko et al. [35], the approximation is also considered in this work for constructing a numerical scheme for solving the fractional differential equation.The numerical approximation for the fractional-order derivative can be derived by employing fractional Taylor series expansion.
This implies

2-Construction of Fractional Numerical Scheme
In the initial level of construction, a difference equation with some unknown parameters is constructed.In a later stage, the unknown parameters will be determined.Also, the scheme will be constructed for the system's first Equation 3. The constructing procedure for Equation 4 will be similar.
The scheme is an explicit and implicit type fractional scheme.The first stage of the scheme is explicit, whereas the second stage is implicit.For its first stage, a difference equation for Equation 3 is expressed as: This stage (Equation 11) is explicit.For the implicit stage, a difference equation in implicit form is constructed as: where , , , , and  are unknown parameters determined later.Using fractional Taylor series for   +1 and      +1 as: The flow chart of the implementation of this scheme is described in Figure 1.

3-Stability Analysis
A Von Neumann stability criterion is employed to check the stability of the constructed scheme.System of Equations 3 and 4, having only diffusion terms, are expressed as: where Applying the first stage of the constructed scheme on Equation 18 using second-order central difference spatial discretization yields: Applying the second stage of the constructed scheme to Equation 18 yields: Substituting Equation 23into Equation 24 is obtained: Equation 25 is simplified to: where  ̅ is an identity matrix of size 2 × 2. Let   be the maximum Eigenvalues of , and then the stability condition is expressed as:

4-Results and Discussions
The proposed scheme for time discretization is third-order accurate in time according to the Taylor series and secondorder accurate in space.Observing the central difference formula's spatial second order can confirm the spatial order.The order of the difference formula can be proved by applying the Taylor series.Also, spatial accuracy can be enhanced by using a higher-order central difference formula for the spatial term.Since the order of the proposed scheme is greater than the first order in space and time, the scheme is consistent according to the considered analysis.Therefore, according to the Lax equivalence theorem, the consistent finite difference scheme will converge if it is stable for linear PDEs.Therefore, the proposed numerical scheme will converge conditionally in the classical case.
The first stage of the proposed scheme is simply the first-order forward Euler scheme that can be capable of finding solutions to fractional non-linear differential equations without linearization.But, the step size should be chosen appropriately so that both stages of the proposed scheme converge.The convergence of the scheme also depends on the chosen values of parameters in the given set of differential equations.Stability is only found for a fractional parabolic equation without having source terms.So, the stability condition does not depend on the source term(s), but it can be extended to fractional time-dependent partial differential equations having source terms.Also, Von Neumann stability can be applied to non-linear PDEs by making them linearized.So, some estimated conditions from linearized PDE can be obtained so that before solving those non-linear equations, the condition on step size can be chosen appropriately.The stability condition also depends on the chosen numerical scheme.The proposed schemes involve more than one parameter, so by choosing the appropriate value of the parameter(s), different numerical schemes can be applied without putting much effort into the computational code.The stability condition may vary for every chosen numerical scheme.Since the proposed scheme is explicit and implicit in two stages, it may require an iterative scheme to solve the difference equation obtained by applying the proposed scheme.So, the solution is found iteratively, and this iterative procedure will stop if the given stopping criterion is met.
Figures 2 to 9 show the results obtained by solving Equations 3 and 4 numerically with some initial and boundary conditions.Figure 2 deliberates the comparison of the proposed scheme with existing numerical schemes.For comparison purposes, one of the existing schemes is backward in time and central in space.This comparison is related to fractional schemes.Since the backward in time and central in space (BTCS) scheme is first order in time, it takes more iterations to converge.The Crank-Nicolson scheme is second-order in space and time, converging faster than BTCS.The proposed scheme is third-order in fractional time, faster than the other two schemes.Figure 3 shows the impact of climate feedback parameter λ on global mean surface temperature  and global mean deep ocean temperature   .The impact of the parameter is shown on both time and space variables.Global mean surface and global mean deep ocean temperatures decrease incrementally in the feedback parameter.Figure 4 shows the impact of deep-ocean heat uptake parameter γ on the global mean surface and ocean temperatures.Both types of temperatures de-escalate by the growth of deep-ocean heat uptake parameters.But the global mean surface temperature has both increasing and decreasing behavior when the heat uptake parameter varies.The impact of heat source parameter Q on the global mean surface and global mean deep ocean temperatures are displayed in Figure 5. Figure 5 shows that global mean surface and global mean deep ocean temperatures escalate by rising heat source parameters.Figures 6 to 9 show the contour plots for global mean surface and deep ocean temperatures.Table 1 compares the results obtained by the proposed and existing schemes with the Matlab solver ode45.The trapezoidal scheme is second-order and does not require any restriction on step size to converge.The Matlab solver ode45 is a ready-made code that can solve any scalar or system of linear and non-linear ODEs.This solver uses a highorder numerical scheme to find the solution(s) to the given ODE(s).The difference of numerical solutions obtained by the proposed/existing schemes and Matlab solver ode45 is found, and then the norm of absolute difference is listed in Table 1.Since the numerical solution is considered in place of the exact solution, the norm of error is still reduced by increasing grid points N for the particular case.

5-Conclusion
This work was related to proposing a numerical scheme for solving fractional and classical time-dependent partial differential equations.We arrived at a more general and realistic conclusion by employing a fractional mathematical model to simulate climate change.The scheme was third-order accurate, and its order can be observed from its construction procedure.The stability analysis of the scheme was also provided.The effect of the contained parameters on the global mean and surface temperature was shown in the graphs.As a result, we concluded that a fractional version of the model produces a satisfactory equivalent result for climate control.From the obtained results, it was observed that increasing feedback parameters de-escalated global mean surface and deep ocean temperatures, and global mean surface and deep ocean temperatures decayed by increasing deep ocean heat uptake parameters.The proposed numerical scheme can be further applied to epidemiological disease models and other problems [40][41][42][43] in fractional calculus.

6-2-Data Availability Statement
The data presented in this study are available on request from the corresponding author.

6-3-Acknowledgements
The authors wish to express their gratitude to Prince Sultan University for facilitating the publication of this article through the Theoretical and Applied Sciences Lab.

6-6-Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this manuscript.In addition, the ethical issues, including plagiarism, informed consent, misconduct, data fabrication and/or falsification, double publication and/or submission, and redundancies have been completely observed by the authors.

Figure 1 .
Figure 1.Flow chart of considered algorithm