{'fraction_non_alphanumeric': 0.040030703386676594, 'fraction_numerical': 0.012606996650539635, 'mean_word_length': 4.324210526315789, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linear Galerkin-Legendre spectral scheme for a degenerate nonlinear and nonlocal parabolic equation arising in climatology

A special place in climatology is taken by the so-called conceptual climate models. These relatively simple sets of differential equations can successfully describe single mechanisms of climate. We focus on one family of such models based on the global energy balance. This gives rise to a degenerate nonlocal parabolic nonlinear partial differential equation for the zonally averaged temperature. We construct a fully discrete numerical method that has an optimal spectral accuracy in space and second order in time. Our scheme is based on the Galerkin formulation of the Legendre basis expansion, which is particularly convenient for this setting. By using extrapolation, the numerical scheme is linear even though the original equation is nonlinear. We also test our theoretical results during various numerical simulations that support the aforementioned accuracy of the scheme.


Introduction
In climate dynamics, one distinguishes between various models according to their complexity and the number of physically resolved phenomena. Overall, there is a hierarchy of models in which on the one end we have General Circulation Models (GCMs), and on the other -Conceptual Climate Models [9,43,45,48]. In between, there are models of Intermediate Complexity that frequently utilize coarser grids or a larger number of parametrizations of various phenomena than GCMs [8]. General Circulation Models are the most advanced descriptions of the global flow of the planetary atmosphere or ocean along with many physical quantities such as pressure, temperature, and water vapour to name only a few. Because they are systems of nonlinear partial differential equations to be solved for long times on a sphere, they usually require supercomputers to analyse and resolve. Thanks to that, they are used to gain quantitative insight about the various interactions between different constituents of the Earth system. Moreover, they are crucial in making predictions of different scenarios concerning, for instance, greenhouse gas emissions and their impact on the mean temperature.
On the other side of the spectrum are Conceptual Climate models that usually are low-order dynamical systems. Their role is to describe several mechanisms of climate along with their nonlinear interdependencies rather than simulate the whole Earth system. They are very useful to investigate the steady states of the climate and their bifurcations which might be difficult to resolve in the full GCMs. For example, conceptual climate models may only focus on energy balance to describe the global or zonally (longitudinally) averaged temperature and steady states of the climate that arise from it. The seminal works concerning this concept have been done by Budyko [4], Sellers [46] for zerodimensional case (global mean), and North [33,34] for meridional evolution of temperature (zonal mean). Overall, in energy balance models, one usually assumes that the temperature distribution T = T (x, t), where x denotes the sine of latitude, is governed by the incoming solar radiation, outgoing infrared radiation, and some mechanism of horizontal (meridional) energy transport. Furthermore, as was first noted in [3] the reflectivity of the surface may depend on the history of the temperature. This leads to a nonlinear and nonlocal parabolic equation that is the main subject of our investigations where the memory operator with a kernel K is defined as Here, d = d(T ) is the possibly nonlinear diffusivity while g is the nonlinear source term. Notice that the nonlocality enters the equation through the latter and the second-order differential operator is degenerate for x = 1. A derivation of the above model in the framework of climate dynamics is given in the next section. In early works, the model was analysed from the climatological point of view in [32] where the natural efficiency of Legendre orthogonal polynomials was noticed. The linear case was solved exactly for several first modes while the nonlinear one, with the diffusivity proportional to the flux as proposed in [49], and iteratively in [25]. In this work the authors noted that nonlinear effects have a significant impact on the sensitivity of the model to variations in the solar constant. Moreover, due to a nonlinear source, which is a consequence of the ice-albedo feedback discussed below, the equation can have several steady states representing different climates. Their stability, bifurcations, and sensitivity to parameter perturbations are of high importance in climatology. The simple energy balance model has been generalized in several ways by adding additional degrees of freedom. For the globally averaged model, one can mention adding a mass balance to account for ice sheet variations [24,20,39,40] or the amount of CO 2 in the atmosphere [18,17]. This is especially relevant for understanding the oscillations of ice ages and their rhythmicity. Further information concerning this topic can be found in [29,9,15,37,41,10]. A modern review of decades of research on energy balance models can be found in a readable book by North and Kim [36]. There is also a broad literature concerning the mathematical aspects of the above problem. Its existence, uniqueness, and regularity of solutions was investigated by Hetzer in [21]. Moreover, an essential problem for applications -parameter estimation -was analysed in [44,5]. As was noted above, steady states of the problem can bifurcate what can lead to hysteresis. Mathematical analysis of this phenomenon was given in [13]. Some other interesting results concerning the local version of the model were obtained by Diaz [11] where the general mathematical theory of energy balance models has been developed for a possibly degenerate nonlinear diffusivity. Furthermore, a generalization of the domain from the 2-sphere into a Riemannian manifold without the boundary was analysed in [2] where also a numerical treatment was conducted. Additionally, since in some parametrizations the source can exhibit a discontinuity in the temperature, a free boundary can arise. We give a climatological background of this phenomenon in the next section and the reader can consult [14] for a thorough treatment. Lastly, we also would like to mention broad studies conducted from the dynamical systems point of view in which the horizontal transport is modelled by a relaxation term for the globally averaged temperature [30,52,53].
There are several accounts of the numerical treatment of energy balance models. The early simulations were based on spectral Legendre decomposition in space and the first-order implicit Euler scheme in time [35]. As was also noted in [33] the model is amenable for such a treatment due to the rapid convergence of the orthogonal series. Recently, several papers concerning the finite element method were published. In [2,1] the situation set on a two-dimensional Riemannian manifold was solved in the case of nonlinear diffusivity modelled by a p-Laplacian. On the other hand, in [22] a finite volume WENO method has been applied to solve the problem where the surface of the Earth is split into land and ocean fractions. Similar analysis but with emphasis on equilibrium solutions was given in [23].
In this paper, we design and analyse a two-step spectral method based on the Galerkin-Legendre approximation that is weighted in time. In particular, our method contains the second order in time Crank-Nicolson scheme. The main motivation of this paper is to present a rigorous convergence treatment of the early ideas of energy balance simulations [35]. Spectral methods are known for their exponential accuracy provided the regularity of the solution. This can significantly reduce the computation time since only a small number of terms in the expansion needs to be calculated. This is especially relevant for dealing with nonlocal operators in time that require all information about the history in each time-stepping iteration. Our approach is based on similar estimates obtained in the finite element setting presented in the classical works of Douglas and Dupont [16]. However, as is seen from (1) our PDE becomes degenerate at x = 1 which causes several difficulties. We overcome them by introducing a weighted L 2 space in which the solution is sought. We also utilize the extrapolation of coefficients in the same way as was done in [28] to obtain a linear method. Galerkin finite element method has also been applied to parabolic equations with nonlocal terms in [6]. On the other hand, Galerkin spectral method has recently been applied to nonlinear parabolic problems in [27]. The reader can find state-of-the-art surveys of all variants of spectral methods in [7,47].
The paper is structured as follows. In the next section, we present a derivation of our main model (1) in the climatological setting. In Section 3, we deal with a semidiscrete scheme where we discretize the spatial variable leaving a continuous dependence on time. The method attains the optimal spectral accuracy. The convergence proofs are given. Furthermore, Section 4 concerns the fully discrete method where a two-step weighted scheme with extrapolated coefficients is constructed and analysed with respect to the convergence. Thanks to the extrapolation, our method is linear despite the fact that the original problem may be highly nonlinear due to the diffusivity. In Section 5, we provide implementation guidelines, solutions that improve calculation performance, and simulations that verify previously proved estimates. We end the paper with a conclusion with some prospects for future work.
In what follows, the generic letter C will be used frequently to denote any positive constant that can depend on the solution and its derivatives but not on the discretization parameters. Moreover, as in the usual practice, the value of C can change even in the same chain of inequalities. This will not cause any confusion or lack of rigour.

Background on climate dynamics
To set the stage for our subsequent reasoning, we review the derivation of the main equations in the setting of climate dynamics. Some more elaborate exposition can be found in [19,40].
We consider a simple conservation of energy written in terms of the zonally and vertically averaged temperature that is symmetric with respect to the equator. Let T = T (x, t) be the mean temperature at time t and latitude θ where x = sin θ. We prefer to use x instead of θ since then, the spherical Laplacian simplifies substantially. The basic model is the following where c is the heat capacity of Earth, R i is the incoming solar radiation that reaches the surface (insolation), R o is the outgoing infrared radiation, and H represents the horizontal transport.
We can further parametrize these various constituents. As the amount of heat Q reaches Earth its fraction α, called albedo, is reflected by the surface. Since, for example, snow (α = 0.80) reflects much more light than the ocean (α = 0.06), we have to distinguish between different types of surface and its albedo. This leads to the so-called ice-albedo feedback that is one of the most important mechanisms regulating the climate. Since ice has a small albedo, it causes more radiation to be reflected, lowering the surface temperature. This, in turn, produces favourable conditions for the formation of ice caps. This phenomenon can be parametrized by letting α to depend on x and T . In the literature one can find many of these functional relations. For example, Budyko proposed that the albedo has two values: one for ice and one for ice-free surface. The boundary between these, called the ice line, depends on the temperature where T i is usually taken as −10 • C. That is, ice starts to form when the mean temperature is below a threshold. This parametrization leads to a free-boundary problem since one has to determine the ice line position x i (t) satisfying T (x i (t), t) = T i . This approach has been analysed, for example in [12].
Other forms of albedo have also been proposed. One of the most common ones assumes continuous dependence on both latitude and the temperature. For example, Sellers suggested a piecewise linear relationship for some choices of T 1,2 . In general, albedo is frequently represented as a bounded, monotone increasing function of the temperature and, possibly, a low-order polynomial in space [24,3,18,40].
We will make such an assumption in the sequel. Moreover, as was suggested in [3] (but also see [42]), the present albedo of the ice-covered ground is determined not only by the actual temperature, but rather by its past values. Whence, we take α to be a function of the nonlocal (memory) operator acting on T , that is where S = S(x, t) is the distribution of insolation over latitude which takes into account an uneven illumination of the spherical Earth. To a good approximation, one can take S(x) where L 2 is the second Legendre polynomial [36]. Furthermore, the absorbed heat is isotropically reradiated into space. This can be parametrized by the Stefan-Boltzmann law (as was done by Sellers) or its linearization R o = A + BT (as was done by Budyko). Since the change in the temperature is relatively small, this simplification is usually justified.
Finally, we have the horizontal transport term that arises due to the uneven temperature distribution along the latitude: heat moves from the equator to the polar regions. To the level of complexity that we want to achieve, we assume that the transport can be modelled by a diffusion term [33] The other approach, by Budyko, is based on using a relaxation term ∝ (T − T ) where T is the global mean temperature. The boundary conditions are taken to model the vanishing energy flux on both the equator x = 0 and poles x = 1, that is we impose −D(u)(1 − x 2 )T x = 0 there. As was found in [32,33] linear diffusion, i.e., D(u) =const., leads to an accurate and robust model and thus, it is an important case to consider. However, as was suggested in [49] and further verified in [25], having the diffusivity being a function of the gradient has a significant effect on the stability of polar ice caps. This parametrization leads to the p-Laplacian operator which, in this case, was analysed in [11]. Putting all obtained terms and renaming the source term, we arrive at (1).
We note that the heat capacity c is assumed to be constant, however, in [44] it was suggested that it may also be a function of the past values of temperature and latitude. Investigating a model of this type along with gradient dependent diffusivity will be the subject of our future work.
3 Spectral discretization with respect to space

Definitions and assumptions
We begin with some preparations. By L 2 (0, 1) we denote the usual Hilbert space of square-integrable functions with a norm · . Similarly, H m (0, 1) where m ≥ 1 is the Sobolev space of m−th weakly differentiable functions with the norm denoted by · m . As for the scalar product, we will only use L 2 one and write (·, ·). We also introduce an intermediate space where the solution to (1) lives Due to the degeneracy of the equation (1) at x = 1, the above weighted L 2 space is a natural choice.
In addition, because of that reason, the corresponding quadratic form needed for the definition of a weak solution is not coercive (it is only weakly coercive). One of the standard ways of dealing with that problem is to introduce a transformation which leads to where f (x, t, u, Ju) = e −t g(x, t, e t u, J(e t u)) and ψ(x, s) = e −t T 0 (x, s). Notice that we should have written D(u) = d(e t u), i.e. indicating the explicit dependence on time. However, since t ∈ [0, t 0 ] we have e t bounded. Therefore, omitting it from the diffusivity will not produce any quantitative effects in the proofs below. We thus commit this slight abuse of notation and abstain from writing explicit time dependence. The reader will see every reasoning below can be repeated for time-dependent cases with essentially no changes. Now, by multiplying the above by v ∈ V and integrating by parts from x = 0 to x = 1 we obtain a weak formulation which is the basis for the Galerkin method where the quadratic form a linear in the second and third argument is defined by We will also write a(u, v) := a(1; u, v) which implies that a(u, u) = u 2 V . Moreover, from now on, if it does not pose any threat to unambiguity, we will suppress writing the independent variables.
Concerning the assumptions imposed on various parameters, we make the natural choices that are also required for existence and uniqueness (see [11,44,21]). Specifically, we assume that the diffusivity D = D(u) and the source f = f (x, t, u, w) are smooth with Moreover, we take the kernel of the nonlocal operator (2) to have a minimal regularity for wellposedness In several places, we will further assume that the solution of (5) is sufficiently smooth, which, in turn, requires more regularity on D and f . We also would like to note that considering a discontinuous source case is one of the subjects of our future work.

Numerical scheme
In the spatial discretization, we use the Galerkin scheme for which the test and trial functions (see [7]) belong to the following finite-dimensional space where P N is the polynomial space of degree N . As for the orthonormal basis where L 2i is the Legendre polynomial of degree 2i. We obviously have (φ i , φ j ) = δ ij . Moreover, this choice is particularly convenient for linear diffusion, i.e., D(u) = const., because it constitutes the eigenbasis for the second-order operator This automatically diagonalizes the stiffness matrix and facilitates computations in the important linear case or makes the matrix sparse for weakly nonlinear diffusion. By weighting (6) with respect to V N we formulate the Legendre-Galerkin numerical scheme. We where ψ N is the appropriate approximation to the initial condition. Note also that we have refrained from writing the x variable as arguments since the inner product is taken with respect to it. This slight abuse of notation should not cause any misunderstandings. Henceforth, we will use the orthogonal which has the following spectral accuracy (for a detailed discussion see [7], formulas (5.4.12) and (5.4.17)) where V (·) denoted the total variation and l ≥ 1. The error of the approximation with P N in V -norm is better than that in the Sobolev space.
with u and u t bounded. Then, for sufficiently large N we have where L is defined in (10).
Proof. We start by writing which follows from the integration by parts. Since Moreover, since φ i is the eigenfunction of L (see (10)) we have However, we can also write φ i = λ −1 i Lφ i which after iteration implies that where the last inequality follows from Plancherel's identity. Further, since , which for sufficiently large N implies the sought error estimate in the V -norm.
We note that the above is not the only optimal choice of a projection that we can use. For the main equation (11) it is much more convenient to use the following elliptic (or Ritz) orthogonal projection which utility in obtaining an optimal order of convergence was discovered in the early days of mathematical finite element analysis [50]. Before we proceed the convergence result for (11) we state several auxiliary lemmas concerning Ritz projection (15). First, we show that it has the optimal order of approximation both in V and in L 2 (0, 1). Lemma 2. Let u(t) ∈ H 2m (0, 1) for m ≥ 1 and each t ∈ [0, t 0 ] with u and u t bounded. Then, for sufficiently large N we have the following estimate in V and L 2 Moreover, if additionally u t ∈ H 2m (0, 1) then the time derivatives of the errors satisfy Proof. By (7) and since R N u ∈ V N we have where v ∈ V N . By the orthogonality of the Ritz projection (15) the last term above vanishes leaving By choosing v = P N u and using Cauchy-Schwarz inequality, we immediately have Furthermore, the application of Lemma 1 leads to which is the first assertion.
To show the L 2 estimate of the error, we follow the duality argument (see [50]). For a fixed u ∈ H 2m (0, 1) let w be the solution of By putting v = w we immediately obtain the stability estimate w V ≤ C u − R N u . Using the definition of a as in (7) and reintegrating by parts, we can obtain where the assumption of bounded D has been used in the first inequality while the stability estimate in the last. Now, by choosing v = u − R N u in (19) and using the definition of the Ritz projection (15) we have Furthermore, we can use Cauchy-Schwarz inequality along with (18) and (14) to infer that The proof of (16) is finished after utilizing (20). The reasoning used in showing the time differentiated version (17) is similar and based on a derivative of the Ritz projection definition (15) a(D(u); From here it follows that

Now, Cauchy-Schwarz inequality yields
. We can make the norm of differences sufficiently small by an orthogonal approximation, i.e., by choosing v = P N u t . Whence, , and by another estimate we can write We can now use the -Cauchy inequality, that is to transform our estimate to Hence, the V -norm part of (17) follows from (14) and (16). The L 2 estimate again follows from the duality argument. Let w be as before in (19) but (u−R N u) t on the right-hand side. Reasoning similarly as before we have with v ∈ V where in the third equality we have moved the derivative according to (21) while in the last we have introduced w. Thanks to that, with v = P N w we can bound the above by previously obtained estimates in the V -norm The first two terms can be tackled exactly in the same way as above, and the third term can be integrated by parts to obtain by stability estimate for elliptic equation (20). Finally, combining the two above estimates with (16) yields (21) and finishes the proof.
As we have seen, the elliptic projection has the same order of accuracy as the standard projection P N . Note also that the error in the V -norm is bounded by N 1−2m for sufficiently regular functions u. This norm involves the first derivative and is similar to the standard H 1 norm with u V ≤ u 1 . This can be compared with the classical result concerning the approximation in H 1 which states that in that case the error is proportional to N 3/2−2m which is larger than the former estimate. This again shows that choosing V as the trial space is very natural and optimal.
The Ritz projection is also bounded in the maximum norm which is shown below.
Lemma 3. Let R N be defined as in (15) and u ∈ H 4 (0, 1). Then, for sufficiently large N we have Proof. We will use the following polynomial inverse inequalities, one for differentiation (Markov inequality, [51], p. 218), and the others for summability (see [7], formulas (5.4.3) and (5. and use the above inverse inequality and the error estimate (16) from which the conclusion follows for sufficiently large N .
The last auxiliary result is the classical Grönwall-Bellman's lemma that we state without the proof (which can be found in [26]).
Lemma 4 (Grönwall-Bellman). Let F (t) and y(t) be continuous, non-decreasing, and nonnegative functions. Then We are ready to prove the main results of this section concerning the semidiscrete numerical scheme. First, we show that the method is stable in time.

Theorem 1 (Stability).
Assume that f = f (x, t, u, w) is bounded with respect to u, and w, that is where g = g(x, t) is smooth. Then, if u N is the solution of (11), we have where D − is given in (8) while λ 1 is defined in (10).
Proof. The proof is standard: in (11) choose v = u N , use Cauchy-Schwarz inequality, and the boundedness assumption to obtain Now, since a(D(u N ); u N , u N ) ≥ D − a(u N , u N ) (see (7) and (8)) we can write which after division by u N (t) leads to The left-hand side is bounded from below by the infimum over the whole space V , hence According to the standard theory of elliptic PDEs, the infimum can be interpreted as where λ 1 − 1 is the smallest eigenvalue of the problem The value of λ 1 is given in (10). Finally, multiplying the inequality by the factor e D − (λ 1 −1)t and integrating yields the sought result.
Again, we can see the close relationship between the linear diffusion generated by the operator L and the fully nonlinear case. We can now proceed to the convergence proof.
Theorem 2 (Convergence). Let u(t) and u N (t) be solutions of (6) and (11), respectively belonging to H 2m (0, 1) with m ≥ 2 for each t ∈ [0, t 0 ]. Moreover, assume that u x and u t are bounded. Then, for sufficiently large N and for each t ∈ [0, t 0 ] we have where the constant C depends on u, u x , u t , D, f , and K.
Proof. We start by decomposing the error Therefore, since we have already proved the estimates on r N in Lemma 2 we have to focus only on e N . To this end, we will write the error equation.
where we have used (11). Furthermore, by using the main equation (6)  Therefore, by taking v = e N ∈ V N we are led to the estimate And from here we are left to estimate the three remainders ρ i . The first one comes from (17), while for the second we use (8) and obtain Now, the gradient of the Ritz projection is bounded according to Lemma 3 and hence The third remainder can be bounded using (8) Now, the nonlocal operator Ju is defined by (2) and hence where in the last inequality we have changed the integration variable s t−s and used the fact that the integral of a positive function over [0, t] is larger than over [t − τ, t]. Finally, we can combine all estimates of ρ i and use -Cauchy inequality where necessary to extract e N 2 V and bound products of norms in terms of the sum of their squares. In effect, we arrive at By integrating the above, we obtain Since K is integrable, its integral is continuous, hence bounded, and and we can invoke Grönwall-Bellman's Lemma (Lemma 4) to arrive at . This, Lemma 2 and the fact that u − u N = r n + e N implies the assertion and completes the proof.

Weighted linear time-stepping scheme
We would like to fully discretize (11) to obtain a time-stepping numerical scheme. However, to reduce the computational cost, we would like to obtain a linear method. To this end, we will use the extrapolation of the nonlinear coefficients (see [50,28]) along with a θ-weighed scheme.
As a preparation, we introduce the temporal grid on [0, t 0 ] with a step h is defined as where N 0 is the number of subintervals. Furthermore, if U n is a grid function, that is a function defined on t n with n = 0, ..., N 0 that can be though as a piecewise constant function on [0, t 0 ], we define the usual finite difference Furthermore, we have to discretize the nonlocal operator (2) and in general the discretization can be written for τ = hM where ρ 0 (h) is the local consistency error, and w i (K) are weights. In particular, we can choose the rectangle rule to have for 0 ≤ i < M or trapezoid rule Provided sufficient smoothness, the orders of the above quadratures are: first, for the rectangle, and second for the trapezoid. Note that we have used the so-called product integration rule, that is, we have discretized the unknown function U leaving the exact integral of the kernel. This procedure guarantees the maximal order of convergence for sufficiently smooth U (see [26]).
To state the method, we define the θ-averaged value U n−θ := θU n−1 + (1 − θ)U n , 0 ≤ θ ≤ 1, and its extrapolation through the past two time steps It can be easily seen by Taylor series expansion that Now, we fix N > 0 as the number of terms in the Galerkin semi-discrete solution of (11) and set U n to be the numerical approximation to u N (t n ). Note that we will omit writing subindex N in our fully discrete approximation. We then propose the following scheme for solving (5) where f h is defined as which is a full discretization of the source nonlinearity. Note that since we are using extrapolation (25) in D and f we are required to solve only a linear system of equations in each time step. This technique of removing nonlinearity is a classical move developed in [28]. Due to the weighted nature of (26) we obtain the first-order backward Euler method for θ = 1 and second order Crank-Nicolson scheme for θ = 1/2. Accordingly, with respect to the requirements, one can choose either the rectangle or a trapezoid quadrature in (22). To complete the numerical scheme, we have to state the initialization procedure. Since the extrapolation produced a two-step method, we have to carefully start the iteration for the order to be preserved. We use the predictor-corrector method. The predictor step is based on setting U 0 = ψ N (0) and solving for W in Then, we correct the value of W by where W 1−θ = (1 − θ)U 0 + θW . This gives us two starting values U 0 , U 1 that can be plugged into the time-stepping (26). We now move to the convergence result.
Theorem 3 (Convergence of the full discrete scheme). Let u(t) ∈ H 2m (0, 1) for each t ∈ [0, t 0 ] with m ≥ 1. Furthermore, assume that u x , u t , and u tt are bounded. Then, where the constant C depends on u, ψ and all their derivatives, and ρ 0 (h) is defined in (22).
Proof. Similarly as in the proof of the semidiscrete scheme, we start by decomposing the error u n − U n = u n − R N u n + R N u n − U n =: r n + e n , where u n = u(t n ) and we keep our convention not to write N explicitly (since it is fixed). Then, we start writing the error equation with v ∈ V N (δe n , v) + a(D(U n−θ );ê n−θ , v) where in the last equality we have used (26). Expanding further, we have where this time we have used the main equation (6). In the above, we would like to put v =ê n−θ in order to derive estimates on the error. Before that, however, note that which can be shown by expanding the definitions of δe n andê n−θ . Note that we see that taking θ = 1/2 kills the O(h) term. Therefore, with the aforementioned choice of v, we can write where the remainders ρ i are understood from (4). We will now bound each of them. We start with the difference in time derivatives δR N u n − u n−θ t ≤ δr n + δu n − u n−θ t . Now, by Taylor expansion at t = t n−θ = (n − θ)h we obtain Whence, by Lemma 2 we have Next, the remainder associated with a can be bounded thanks to (8) and the difference between extrapolated and the value of the exact solution can be estimated as follows since, by construction, u n−θ approximates u n−θ to the second order. Next, we proceed with the third remainder to obtain where we have used the -Cauchy inequality to extract the V -norm of the error. However, exactly as above, due to Taylor expansion at t = t n−θ = (n − θ)h we have for any function y of time To apply this estimate to the ρ 3 we have to show that the gradient of the second time derivative of Ritz projection is bounded in V . To this end, differentiate the definition (15) twice to obtain for sufficiently large N by Lemma 2. Whence, according to Lemma 2 we finally obtain The last remainder involves a nonlocal operator. First, we estimate the difference by the regularity assumption on f

Now, by (22) we have
by integrability of the kernel K. Hence, Putting all estimates of ρ i that we have obtained so far and defining ρ(N, h) := ρ 0 (h) + h 2 + N −2m , leads to where we have again used the -Cauchy inequality to extract the L 2 norm of the errors. Now, since K is integrable, the weights w i (K) are bounded which along with a simple real number inequality (a + b) 2 ≤ 2(a 2 + b 2 ) yields a nonlocal recurrence By iteration, we then have and therefore, Now, we are left with estimating the error made in the first predictor-corrector step (28)-(29) since the sum above involves the initial condition. To this end, set g 1 = W − R N u 1 . Similarly as above, we can show that for the predictor (28) we have Then, U 0 − u 1−θ 2 ≤ e 0 2 + r 0 2 + Ch ≤ C e 0 2 + N −2m + h , and it follows that Next, we move to the corrector stage (29) to obtain δ e 1 2 ≤ C W 1−θ − u 1−θ 2 + ρ(N, h) 2 .
Reasoning as above, we have Using the estimate on g 1 from the prediction stage, we obtain Now, by going back to the estimate on the finite difference of the error, we finally have The estimate of the initial value errors e −i where i = 0, 1, 2, ..., M follows from what ends the proof.

Implementation and numerical illustration
Now we are concerned about the practical use of the aforementioned algorithm and its efficient implementation. Let U n be expanded into our basis (9) then, plugging into (26), by orthonormality of {φ i } i we obtain where y n = {y n i } N i=0 is a vector of solutions, while the stiffness matrix A = {A ij } N i,j=0 and the load vector In the linear case, the stiffness matrix is diagonal A = diag(λ i ) i with eigenvalues (10). The implementation requires solving (30) in each time step for y n which reduces to inverting the nonsingular matrix I + (1 − θ)hA(U n−θ ). For the linear diffusion, this matrix is constant over time and the inversion needs to be done only at the initialization phase. We would like to numerically verify the above results concerning convergence. To this end, we will calculate the order of the method for the Crank-Nicolson scheme with θ = 1/2, with trapezoidal quadrature (24) and two choices of nonlocal operator kernels. The first one is the Gaussian as was suggested in the original work [3] where A is the amplitude and σ 2 is the variance. This kernel is responsible for a short memory effect due to the exponential decay of its tail. In contrast to that, we can also consider a long memory kernel obtained by the power law describing the heavy tail In the above, the prefactor involving the gamma function has been chosen to be consistent with the Liouville fractional integral which arise after substitution s → t − s and allowing for the infinite memory with τ → ∞. Fractional integrals and fractional derivatives are important in many applications and there is a very vigorous research going on this topic: both applied and pure. The interested reader can find additional information in [31].
In our numerical examples, we would like to test two cases: either linear diffusion with nonlocal operator or nonlinear diffusion without nonlocality. Errors are calculated in the L 2 norm and the coding is done in Julia programming language. It is a free high-performance and high-level dynamic programming language that, apart from many others uses, performs very well in numerical simulations where one has to conduct a large scale computations. Julia neatly combines ease of coding with speed of execution (for an introduction for scientists see [38]). We have implemented several performance mechanisms: • All x-integrals are calculated using Gaussian quadrature with pregenerated Legendre weights.
• Where possible, we utilize parallel computing with several threads/cores (multi-threading). For example, the stiffness matrix and load vector entries, as well as quadrature weights (23)- (24) can readily be calculated effectively in this way.
• Quadrature weights need only be generated once as soon as the kernel K and the time step h are fixed.
• For linear diffusivity, the stiffness matrix is pregenerated. For nonlinear diffusion, we utilize multi-threading.
We have found that the approach based on parallelism is highly efficient with respect to the single core calculations. Below, for concreteness we set τ = 0.4 and t 0 = 0.5. Moreover, the initial condition is always taken to be ψ(x, s) = cos(πx) 1 + s , which models a high temperature at the equator and low at the pole. Overall, we solve three cases of our problem which test the effects of nonlinearity and nonlocality in several ways. The first example introduces pure nonlinear effects into the equation: both in the diffusivity and the source. The former has been chosen to satisfy positive definiteness while the latter to ensure three nontrivial steady states that are usually present in energy balance models. They describe (stable) cold, (unstable) intermediate, Figure 1: Numerically calculated L 2 error between solutions with N max = 30 for different N and problems given in (32).
and (stable) warm climates (more information can be found for ex. in [19]). Next, two examples investigate the notion of history: short (Gaussian) and long memory (fractional integral) kernels.
For the latter, we choose a representative value α = 0.5 since the results for others do not present any significant qualitative differences. First, to check the spectral accuracy of the spatial discretization, we fix the time step to be h = 2 × 10 −2 and compare solutions obtained for different N with the one calculated for N max = 30. The latter is treated as "exact". The time step h can also be chosen smaller and in that case the error would decay even faster. The comparison is done at t = 0.125 to not let diffusive effects to force the solution to decay, which can happen for larger times. Results of our simulations are shown in Fig. 1 where a semi-log plot of errors in presented. As we can see, the spectral accuracy is confirmed for all considered cases, that is the logarithm of the error behaves approximately linearly. The error saturates at machine epsilon (that is M = 2.2 × 10 −16 ) for N ≈ 11 − 13 for all problems which means that using several degrees of freedom can grant the (numerically) perfect accuracy. We can also see that the convergence to zero is slightly faster for no or short memory examples. The error could be made even smaller if we decreased the time step, however, we wanted to clearly show the error trend.
The temporal order is calculated in a similar way. Here, we use N = N max = 20. As we have seen in above computations, this number of degrees of freedom is sufficient to grant a complete spatial accuracy. Our "exact" solution is then precalculated for h min = 2 × 10 −3 and compared with calculations done with larger steps. The time of comparison is taken to be t 0 . In Fig. 2 we present the log-log plot of the L 2 error for different values of h. As we can see, for all considered cases the order of convergence is confirmed to be equal to 2, that is the lines are parallel to h 2 for large values of h −1 (small h). Note that the weakly singular kernel in the fractional integral (31) does not impair the accuracy. This is because we have chosen to integrate exactly the kernel using the product integration rule. The only possible slight loss of accuracy is seen in the nonlinear diffusivity example. However, the graph still confirms the method's second order.
Lastly we evaluate the relative speed of our numerical method. This benchmark is based on comparison with the standard Finite Element Method with second-order interpolation in space and the method of lines adaptive implicit Runge-Kutta time integration. To provide the state-of-the-art computations, we have used the Wolfram Mathematica symbolic computation environment. The specific numerical routine was then NDSolve with appropriate options. To provide some objective, i.e., CPU independent data, we compute the ratio of the computation time for the FEM method and our spectral scheme required to obtain a fixed level of absolute error. That is, for a given tolerance, we calculate the following quantity ratio = time for FEM time for spectral method (33) by choosing the minimal spatial degrees of freedom to obtain a given error. The problem that we solve is the following D(u) = 1 + u, f (x, t, u, w) = u(1 + u), u(x, 0) = cos(πx).
Our results are presented in Tab. 1. The superiority of using the spectral method is evident. As can be seen, it is 40 − 100 times faster than the FEM of choice. This advantage comes from the fact that spectral methods are best suited for finding smooth solutions. Note that, if we were to expect less regular solutions, the FEM would be the method of choice.

Conclusion
Motivated by the efficient use of the spectral method for solving equations of energy balance models, we have provided convergence proofs for Galerkin-Legendre scheme. Spectral approximation to the solution of the diffusive energy balance model has been known since its early days. However, throughout many years despite the successful use of the Legendre method, there has not been any rigorous analysis of convergence and stability. The model brings out a difficulty originating in degenerate diffusivity. Our treatment introduced some special weighted Sobolev space that made the convergence analysis feasible along with optimal estimates on the error. The use of spectral method is also beneficial when it comes to the treatment of nonlocal (memory) operators. Their numerical evaluation requires substantial computing power due to the history of the process. Since spectral methods have exponential accuracy it is then possible to reduce the number of degrees of freedom on the spatial part of the problem and use the remaining computational power to treat the memory effects on the temporal side. Discretization in time using extrapolated coefficients grants a linear second-order method that can quickly compute the solution to the desired accuracy. Moreover, using a parallelized code fragment helped to improve the performance even further. Due to the nonlocal operator that requires to take into account the history of the evolution, the initial computational cost of the simulations can be high. Thanks to the spectral accuracy and multithread computations, we were able to reduce it. In future work, we will extend our method to some more general equations with nonlocal specific heat and diffusivity proportional to the gradient. That is to say, we plan to consider the following problem suggested, for example, in [11] c(x, T, JT )T t = (d(T )(1 − x 2 )|T x | p−2 T x ) x + g(x, t, T, JT ), 1 < p < ∞, x ∈ (0, 1), with Neumann boundary condition and the usual choice of the initial condition and the memory operator defined in (2). Two difficulties arise here: one due to nonlinear and nonlocal heat capacity, and the other due to doubly degenerate diffusivity. The latter means that apart from degeneracy at x = 1 we also have to carefully deal with the situation when T x ≈ 0. Both of these generalizations may have a significant impact on the numerical analysis of the problem.
It will also be both relevant and interesting to consider the free-boundary problem frequently found in conceptual climate models [12]. As we have mentioned in Section 2, following Budyko, one would like to introduce the concept of the ice line x i (t) separating the regions neighbouring with T = T i , where T i is a prescribed temperature. The archetypal albedo could then be (3). An additional Stefan-like condition would then be the continuity of the flux at the free boundary x = x i (t). In that case, our spectral scheme could then be coupled with numerical approaches to solving free-boundary problems. Investigating both of these situations will contribute to our programme of a rigorous analysis of climatically relevant nonlinear and nonlocal models.