Next Article in Journal
Node and Network Entropy—A Novel Mathematical Model for Pattern Analysis of Team Sports Behavior
Next Article in Special Issue
On the Existence of Solutions of a Two-Layer Green Roof Mathematical Model
Previous Article in Journal
Finite-Time Attitude Fault Tolerant Control of Quadcopter System via Neural Networks
Previous Article in Special Issue
Mathematical Model of Decomposition of Methane Hydrate during the Injection of Liquid Carbon Dioxide into a Reservoir Saturated with Methane and Its Hydrate
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Approach of the Equilibrium Solutions of a Global Climate Model

1
Departamento de Ingeniería Geológica y Minera, ETS de Ingenieros de Minas y Energía, Center for Computational Simulation, Universidad Politécnica de Madrid, Calle Ríos Rosas, 28003 Madrid, Spain
2
Departamento de Matemática Aplicada, ETS de Arquitectura, Center for Computational Simulation, Universidad Politécnica de Madrid, Av. Juan de Herrera, 28040 Madrid, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(9), 1542; https://doi.org/10.3390/math8091542
Submission received: 21 August 2020 / Revised: 6 September 2020 / Accepted: 7 September 2020 / Published: 9 September 2020
(This article belongs to the Special Issue Modeling and Numerical Analysis of Energy and Environment)

Abstract

:
We consider a coupled model surface-deep ocean effect, where an Energy Balance Model (EBM) is used for modelling the surface temperature and a two-dimensional heat equation represents the evolution of the temperature of the deep ocean. Although the model under study is based on that proposed by Watts & Morantine (1990), here we consider a modified model that incorporates other processes, such as the nonlinear diffusion and the action of coalbedo, depending on the temperature. The stationary states of the model under study, taking the solar constant as the parameter, are numerically attained. The results of the simulation are depicted in a { ( Q , u ) } plot where u is the temperature in the surface and Q is the solar constant. The numerical solution is achieved by means of a finite volume scheme with Weighted Essentially Non-Oscillatory (WENO) reconstruction in space and third order Runge-Kutta scheme, which verifies the Total Variation Diminishing (TVD) property, for time integration. The equilibrium states are accomplished by evolving in time the numerical solution until the stationary solutions are reached. The main novel results of this work concern the numerical obtention of the stationary solutions of both the EBM and the coupled model EBM-deep ocean and the agreement of these results with the theoretically obtained in previous works, where an interval of values of the solar constant Q was obtained with the existence of at least three stationary solutions. In this work, we have numerically obtained more than three stationary solutions for such interval of Q.

Graphical Abstract

1. Introduction

Climate system is very complex and it involves many components and complicate mechanisms. The full behaviour of the climate is something still far to be understood. The pioneering climate Energy Balance Models (EBMs) were introduced by Budyko and Sellers separately (1969) in [1,2], respectively. EBMs are diagnostic models that describe the evolution of the climate for relatively long scales. Several aspects of the mathematical treatment of different versions of climate EBMs have been studied by many authors, such as [3,4,5,6,7], just to name a few of them. In this work, we are concerned with a two-dimensional climate model (latitude-depth), which represents the coupling: mean surface temperature with ocean temperature. It is based on a climate model proposed by Watts and Morantine [8], which is composed of a parabolic equation in a global ocean with a dynamic and diffusive boundary condition, but incorporating additional terms, such as nonlinear diffusion and the coalbedo, depending on the temperature (see [9,10]). Regarding the numerical simulation of this type of climate models, several references can be mentioned, such as [11,12] using finite element methods or [13,14] based on finite volume methods (FVMs). Finite volume methods are very efficient techniques for the numerical approach of the solution of problems governed by balance laws, as those considered in this work. In FVMs, the REA (Reconstruction-Evolution-Averaging) algorithm, see [15], is applied. Thus, this numerical technique involves the computation of integral averages of the solution for each grid cell (also denominated as cell averages) at each time step. Certain reconstruction process is required in order to compute values and gradients from the cell averages where needed. This reconstruction function is usually piecewise polynomial of certain degree, and must be conservative in each one of the control volumes that conform the stencil, so as to keep the conservation properties of the finite volume method. We note that, with the purpose of developing numerical schemes of order greater than one, it is necessary to implement a nonlinear numerical approach, which allows to overcome Godunov’s theorem [16], which states that monotone behavior of a numerical solution cannot be assured for linear finite-difference methods with more than first-order accuracy. As a first approach, second-order accurate schemes can be obtained based on nonlinear piecewise linear reconstruction, such as the minmod method, the Nessyahu-Tadmor scheme, WAF (Weighted Average Flux) scheme, or the MUSCL-Hancock approach (which is second order in space and time), just to mention some possibilities. Complete and detailed descriptions on finite volume techniques can be consulted in [15,17,18]. The first successful attempt to obtain high order monotone numerical schemes are the Essentially Non-Oscillatory (ENO) methods, which were put forward in the pioneering work of Harten, Engquist, Osher, and Chakravarthy in 1987 [19], where the stencil of the reconstruction polynomial is chosen in such a way that it avoids non-physical oscillations in the numerical solution. It is a nonlinear scheme in the sense that it uses known values of the solution in order to select the stencil. Another alternative, widely used nowadays, are the so called Weighted Essentially Non-Oscillatory (WENO) schemes introduced in [20] which, instead of choosing the "smoothest" stencil, as in the ENO approach, a convex combination of all candidate stencils is performed, considering the values taken by certain smoothness indicators, in order to achieve the essentially non-oscillatory property. This WENO procedure is the one implemented in this work to accomplish the high-order nonlinear reconstruction (see, e.g., [21,22,23,24,25] for a detailed description on this technique). Once the reconstruction function is obtained, values of the solution and gradients are attained where needed. In a general case, if we wish to achieve an order of accuracy ( 2 r 1 ) , then we must consider r candidate stencils, each one of them with r cells. A variant of classical WENO approach is the so-called Central WENO (CWENO) scheme [26,27,28,29,30], which makes use of polynomials of different degrees. Other relevant references with different ways to implement WENO methodology can be found in [31,32,33]. With regard to the evolution stage, we remark that we are producing a semi-discrete approach, in which the spatial discretization, performed by means of the FVM, converts the original partial differential equation (PDE) into an ordinary differential equation (ODE). The numerical solution of the ODE requires to implement an ODEs solver. In this work, we have resorted to the third order Runge-Kutta Total Variation Diminishing (RK3-TVD) scheme. Pioneering works on Runge-Kutta TVD methods with different orders of accuracy can be found in [34,35]. The RK3-TVD is a one-step method developed in three stages verifying the remarkable TVD property, which allows for avoiding non-physical oscillations in the numerical solution and it was introduced in the context of second order finite difference schemes by Ami Harten [36]. This property means that being u i n the averaged numerical solution for cell i and time t n , it is verified that T V ( u i n ) T V ( u i n + 1 ) , where T V ( · ) represents the so-called Total Variation, which is given by T V ( u ) = i ( u i + 1 u i ) . The aim of this work is to obtain numerically, by means of the finite volume method with WENO reconstruction and RK3-TVD scheme for time integration, the equilibrium states and an approximation of a subset of the bifurcation diagram in a global climate model taking the solar constant as the parameter. We are interested in computing the stationary states of the EBM, without influence of the deep ocean, and also of the coupling EBM-deep ocean. The reason behind the consideration of these two situations is that both cases are relevant in climatic modelling. In the case of EBMs, without ocean effect, they are frequently used in many studies of global climate, see, for instance [3,4,7,37,38] and references therein, as they allow to predict the average surface temperature of the Earth due to incoming solar radiation, emission of radiation, energy absorption, and greenhouse effects, among other phenomena. On the other hand, the inclusion of the ocean effect has an important influence on the temperatures distribution due to the well-known thermostatic effect of the ocean, absorbing solar radiation and releasing it. Motivated by the importance of both models, in this work we start by dealing with Energy Balance Models (EBMs), in which the evolution of the temperature is governed by a nonlinear diffusive PDE, with the purpose of accomplishing the stationary solutions by evolving in time the numerical solution of the evolution model. Next, we will obtain an approximation of the stationary solutions and part of the bifurcation diagram of a coupled model surface/deep ocean, where an EBM is used to model the surface, and a 2D heat equation is used for the deep ocean. A comparison of the approximated bifurcation diagrams in both cases is also performed. Previous relevant works in which the number of stationary solutions was studied are [39,40,41,42,43]. In the work [44], it is studied the sensitivity of a climatology model with respect to variations in the solar constant and showed the existence of a continuous and unbounded S-shaped set { ( Q , u ) } where Q is the solar constant and u is the temperature. In the present work the effect of the ocean is also considered, which is a novel approach with respect to the aforementioned works. We show numerical results in { ( Q , u ) } plots of the bifurcation diagrams obtained and some of the stationary solutions achieved. The theoretical results stated and proved on multiplicity of steady states in the aforementioned references are numerically verified and completed in the present work and we obtain new information about an approximation of a subset of the bifurcation diagram. The rest of the document is organized, as follows. In the next section, we give a detailed description of the Energy Balance Model considered, and we then introduce the coupled model EBM-deep ocean. The following section is devoted to a brief description of the numerical approach applied, based on the finite volume method with third order Runge-Kutta TVD with WENO reconstruction. The next part deals with the numerical results. Finally, some conclusions and discussion are given.

2. The Energy Balance Model (EBM)

From the mathematical point of view, the 2D EBM (latitude-longitude) has a spatial domain given by a Riemannian manifold M without boundary, simulating the Earth surface, see for instance [9]
C ( x ) u t d i v ( k ( x ) | u | p 2 u ) + R e ( x , u ) = R a ( x , u ) i n M × ( 0 , T ) , u ( x , 0 ) = u 0 ( x ) i n M ,
where C ( x ) is the heat capacity, k ( x ) is the thermal conductivity, R e is the emitted energy, R a represents the absorbed energy, and the exponent p is a real number that is usually taken as p 2 (the case p = 2 is the linear one), u is the gradient of the temperature u. Usually it is taken p = 3 ([45]). In this equation the subscript t denotes partial time derivative. We shall use this notation, very often, hereafter throughout the text. In one-dimensional models the temperature is assumed to be constant over each latitude, getting the problem
C ( x ) u t ( k ( x ) ( 1 x 2 ) p / 2 | u x | p 2 u x ) x + R e ( x , u ) = R a ( x , u ) i n ( 1 , 1 ) × ( 0 , T ) , ( 1 x 2 ) p / 2 | u x | p 2 u x = 0 , x 1 , 1 , u ( x , 0 ) = u 0 ( x ) , x ( 1 , 1 ) ,
where x is the sine of the latitude and we have introduced the weight ( 1 x 2 ) p / 2 after changing to spherical coordinates in (1). Concerning the absorbed energy, R a , this can be obtained as
R a = Q S ( x ) β ( u ) ,
where β ( u ) is the planetary coalbedo, representing the fraction of radiation which is absorbed. It is dependent on the temperature. The coalbedo changes around a characteristic temperature,  10 °C. It is assumed to be discontinuous in the Budyko’s model and continuous in the Sellers’ model. We note that, in the PDE, the nonlinearity of this type (Budyko’s type) is given by a multivalued graph (Heaviside’s type). Concerning the solar constant (Q), it is a positive number that represents the amount of energy received by the Earth per unit of time and unit of space. Although it is not constant, it is considered to be approximately 1360 W/m2. Because the surface of the Earth is four times larger than the cross sectional area, the amount of energy that is distributed throughout the surface is 340 W/m2. Finally, S ( x ) is an insolation function, which allows for the distribution of solar radiation as a function of latitude. It is usually given by polynomial functions. In this work, we have used the expression S ( x ) = ( 5 x 2 ) / 4 to represent less radiation in the Poles and higher radiation in the Equator. As for the emitted energy, we can choose between Sellers’ model, with R e ( u ) = σ u 4 , which utilises the Stephan–Boltzmann’s law, or the Budyko’s model R e ( u ) = B u + A , which makes use of the Newton’s cooling law.
The mathematical analysis of (2) can be found in [46]. The existence of solutions of the problem (1) under the previous hypotheses is proved in [9] by fixed point arguments. One of the model’s main features is its high sensitivity front variation of parameters. The multiplicity of steady states depending on the parameter Q was studied in [40]. There is a bounded interval of Q such that for every Q in that interval, problem (1) has at least three stationary solutions. In [44], the existence of a S-shaped bifurcation branch for the associated stationary problem was proved. The stabilization of the evolution solution of (1) when t + , to a solution of the associated stationary problem of (1) was proved in [40]. Many works are devoted to the mathematical treatment of global climate energy balance models (one layer), among them, we mention [3] and the references therein, and also [4,5]. In [11], a finite element approach is given to a two-dimensional (2D) climate energy balance model without deep ocean effect.

3. The Coupled Model: Deep Ocean-EBM

The mathematical model that we are dealing with is based on the proposed by Watts and Morantine [8] and completed in [13]. This model represents the evolution of temperature within an ocean of depth H. The spatial variables ( x , z ) are x = s i n ( l a t i t u d e ) and z = d e p t h . The spatial domain considered is Ω = 1 , 1 × 0 , H with boundary Γ = Γ H Γ 0 Γ 1 Γ 1 where
Γ H = ( x , z ) Ω ¯ : z = H , Γ 0 = ( x , z ) Ω ¯ : z = 0 Γ 1 = ( x , z ) Ω ¯ : x = 1 , Γ 1 = ( x , z ) Ω ¯ : x = 1 .
One simplifying hypothesis of this model is that a constant temperature on each parallel is assumed. The equation representing the evolution of the temperature in the deep ocean is
U t ( K H R 2 ( 1 x 2 ) U x ) x K V U z z + ω U z = 0 i n Ω × ( 0 , T ) ,
where U ( x , z , t ) is the temperature within the ocean, ω is the vertical velocity, K V is the vertical diffusivity, K H is the horizontal diffusivity, and R represents the radius of the Earth. We assume, as boundary conditions, the following expression for the ocean bottom
ω x U x + K v U z = 0 , o n Γ H × ( 0 , T )
and for the upper boundary (surface of the ocean), we have
D u t D K H 0 R 2 ( 1 x 2 ) p 2 u x p 2 u x x + B u + A + K V U n + ω x u x = 1 ρ c Q S ( x ) β ( u ) o n Γ 0 × ( 0 , T ) ,
where u ( x , t ) is the temperature on the upper boundary, ω is the velocity, K V is the vertical diffusivity, K H 0 is the horizontal diffusivity, R is the radius of the Earth, D is the thickness of the mixed layer (the ocean layer right below the surface, where the relevant heat exchange processes between atmosphere and ocean take place), ρ is the density, c represents the specific heat coefficient, β ( u ) is the temperature dependent coalbedo, Q is the solar constant, B u + A is the cooling term according to Newton’s law, and S ( x ) is the insolation. The full final system representing the coupled model: energy balance on the surface-deep ocean reads
U t ( K H R 2 ( 1 x 2 ) U x ) x K V U z z + ω U z = 0 in Ω × ( 0 , T ) , ω x U x + K V U z = 0 in Γ H × ( 0 , T ) , D u t D K H 0 R 2 ( 1 x 2 ) 3 2 u x u x x + K V U n + ω x u x + B u + A 1 ρ c Q S ( x ) β ( u ) o n Γ 0 × ( 0 , T ) , U x = 0 , on Γ 1 × 0 , T Γ 1 × 0 , T , U ( x , z , 0 ) = U 0 ( x , z ) , i n Ω , u ( x , 0 ) = u 0 ( x ) , i n Γ 0 ,
where we have added initial conditions and homogeneous Neumann boundary conditions at the boundaries Γ 1 and Γ 1 . We note that U ( x , 0 , t ) = u ( x , t ) and U | Γ 0 = u . We have also introduced p = 3 as proposed by Stone [45]. In previous works, [13,14,47], the numerical approximation of the coupled model has been carried out. Furthermore, in [13], latent heat and delay effects are also considered, while, in [47], a land-sea distribution is incorporated to the model. The following hypotheses hold (see [13] for details)
(H1)
β is a bounded maximal monotone graph of R 2 , such that there exist two constants 0 < m < M and ϵ > 0 , such that β ( r ) = m if r < 10 ϵ , β ( r ) = M if r > 10 + ϵ .
(H2)
S L ( Γ 0 ) and S 1 S ( x ) S 0 > 0 a.e. x Γ 0 .
(H3)
w C 1 ( Ω ) .
(H4)
The constants B, A, R, Q, K V , K H , and K H 0 are positive.
We note here that the graph
β ( r ) = m r < 10 , m , M r = 10 , M r > 10 ,
verifies the hypothesis (H1).

Stationary Model

We consider the stationary problem that is associated to (8), which reads
( K H R 2 ( 1 x 2 ) U x ) x K V U z z + ω U z = 0 i n Ω , ω x U x + K V U z = 0 i n Γ H , D K H 0 R 2 ( 1 x 2 ) 3 2 u x u x x + K V U n + ω x u x + B u + A 1 ρ c Q S ( x ) β ( u ) o n Γ 0 , U x = 0 , on Γ 1 Γ 1 .
Concerning the multiplicity of stationary solutions of (10) (depending on Q), see Díaz, Tello [39], where under the hypotheses ( H 1 ) ( H 4 ) , for β defined in (9) and the technical assumption:
10 B + A > 0 w i t h A > 0 , B > 0 , M m S 1 S 0 ,
this result was proved:
(i)
i f 0 < Q < 10 B + A M S 1 ρ c , the problem (10) has a unique solution,
(ii)
i f 10 B + A M S 0 ρ c < Q < 10 B + A m S 1 ρ c , the problem (10) has at least three solutions,
(iii)
i f 10 B + A m S 0 ρ c < Q , the problem (10) has a unique solution.
The fact that the EBM is coupled with the deep ocean model as a diffusive boundary condition and, since the existence of a S-shaped bifurcation branch for the EBM was proved in [44], we expect that ‘a kind of’ S-shaped bifurcation could be obtained for the coupled stationary problem (10). To our knowledge, the bifurcation diagram for the coupled model has not been studied so far. In this paper we obtain some information about the { ( Q , u ) } for the different values of Q being u = U | Γ 0 with U an approximated solution of the coupled stationary problem (10). Concerning the number of steady states, we have obtained more than three numerically approximated steady states for a range of values of Q, where the topological methods proved the existence of at least three (see [39]). We proceed solving the evolution problems until the stationary solution is attained in order to numerically obtain the steady state solutions for the different values of Q. The numerical approach followed in this work is based on a finite volume scheme with third order Runge-Kutta TVD scheme for time integration and using a WENO technique for spatial reconstruction.

4. Numerical Resolution

It is useful to write the problem as a balance law in order to apply the finite volume method. Thus, we rewrite the EBM Equation (7), as
u t f x , u ( x , t ) , u x ( x , t ) x = σ ( x , u ( x , t ) , U n ( x , 0 , t ) ) ,
with the flux
f x , u ( x , t ) , u x ( x , t ) : = K H 0 R 2 ( 1 x 2 ) 3 / 2 u x ( x , t ) u x ( x , t ) w D x u ( x , t ) + K V U z ,
and the source term given by
σ ( x , u , U n ) : = 1 D ( A + Q ρ c S ( x ) β ( u ) + ( ω + x ω x B ) u ( x , t ) K V U n ) .
Concerning the deep ocean, Equation (5), we can introduce the formulation
U ( x , z , t ) t ( F ( x , U x ( x , z , t ) ) ) x ( G ( U ( x , z , t ) , U z ( x , z , t ) ) ) z = Ξ ( x , U ( x , z , t ) ) ,
where
F x , U x ( x , z , t ) : = K H R 2 1 x 2 U x ( x , z , t ) ,
and
G ( U ( x , z , t ) , U z ( x , z , t ) ) : = K V U z ( x , z , t ) w U ( x , z , t ) ,
are the fluxes in horizontal and vertical directions respectively. Additionally, we have the source term
Ξ ( x , U ( x , z , t ) ) : = ω z U ( x , z , t ) .
We solve the problem by means of a numerical scheme that was built under the finite volume framework, with a seventh-order dimension-by-dimension WENO approach and a third order Runge-Kutta TVD scheme for time integration. As aforementioned, the WENO reconstruction applied here makes use of entire polynomials, as introduced, for instance, in [24,25], instead of the more classical pointwise WENO reconstruction ([20,48]). This is especially relevant when solving reaction-diffusion problems where gradients of the solution are involved. The general procedure at each time step consists of solving the Equation (7) getting the cell averages u i n + 1 and use these values as Dirichlet boundary condition to be used in (5) in order to obtain the cell averages U i n + 1 .
We integrate the Equation (7) over the space-time control volumes I i = [ x i 1 2 , x i + 1 2 ] , dividing by the length to get
d u i ( t ) d t = 1 Δ x i f i + 1 2 f i 1 2 + σ i ( t ) l i ( u ( t ) ) ,
where
u i ( t ) = 1 Δ x i x i 1 2 x i + 1 2 u x , t d x ,
is the cell average of u ( x , t ) over the control volume I i and
f i + 1 2 = f u x i 1 2 , t , u x x i + 1 2 , t ,
is the right intercell numerical flux, while
σ i ( t ) = 1 Δ x i x i 1 2 x i + 1 2 σ ( x , u , U n ) d x ,
is the integral average of the source term.
We focus now on the Deep Ocean Model (DOM), which is integrated over the control volume I i j = [ x i 1 2 , x i + 1 2 ] × [ z j 1 2 , z j + 1 2 ] to yield
d U i , j d t = 1 Δ x i F i + 1 2 , j F i 1 2 , j + 1 Δ z j G i , j + 1 2 G i , j 1 2 + Γ i , j L i , j ,
where
U i , j = 1 Δ x i Δ z j z j 1 2 z j + 1 2 x i 1 2 x i + 1 2 U ( x , z , t ) d x d z ,
is the cell average of the temperature of the DOM, whereas
F i + 1 2 , j = 1 Δ z j z j 1 2 z j + 1 2 F x i + 1 2 , U x x i + 1 2 , z , t d z , G i , j + 1 2 = 1 Δ x i x i 1 2 x i + 1 2 G U x , z j + 1 2 , t , U z x , z j + 1 2 , t d x ,
are the spatial integral average of intercell fluxes and
Γ i , j ( t ) = 1 Δ x i Δ z j z j 1 2 z j + 1 2 x i 1 2 x i + 1 2 Ξ ( x , U ( x , z , t ) ) d x d z ,
is integral average of the source term.

4.1. Spatial WENO Reconstruction

We make use of Weighted Essentially Non Oscillatory (WENO) reconstruction from cell averages in order to obtain values of the solution and derivatives where they are needed. In the 2D case, we apply a dimension-by-dimension reconstruction process, as reported for instance in [23,24,25], just to name a few of them. This way to proceed is computationally less costly than the fully 2D WENO reconstruction. There are several variants of WENO approach, among them we can mention the Central WENO scheme, see e.g., [29,49], which considers polynomials of different degree, or the recently introduced WENO scheme with Unconditionally Optimal Accuracy, see [32]. Below, we briefly describe the process followed in this work, using reconstruction polynomials of a general degree M . The description refers to the dimension-by-dimension 2D case.
For each particular cell I i , j , a set of 1D stencils is considered, which, for each Cartesian direction, are given by
S i , j s , x = e = i L i + R I e , j , S i , j s , y = e = j L j + R I i , e ,
where L and R represent the spatial extension of the stencil to the left and to right, respectively. According to the methodology described in [24,25], in the case of odd order schemes, three candidate stencils are taken into account while, for even order, four candidate stencils are considered. Thus, in the case of odd order schemes (that is, even polynomial degrees M), we have
  • One central stencil: s = 1 , L = R = M / 2
  • One fully biased to the left stencil: s = 2 , L = M , R = 0 ,
  • One fully biased to the right stencil: s = 3 , L = 0 , R = M ,
whereas, for even order schemes (that is, odd polynomial degrees M)
  • Two central stencils: s = 0 , L = f l o o r ( M / 2 ) + 1 , R = f l o o r ( M / 2 ) and s = 1 , L = f l o o r ( M / 2 ) , R = f l o o r ( M / 2 ) + 1
  • One left-sided stencil: s = 2 , L = M , R = 0 ,
  • One right-sided stencil: s = 3 , L = 0 , R = M .
The following coordinates transformation is introduced [ x i , x i + 1 ] × [ y j , y j + 1 ] [ 0 , 1 ] × [ 0 , 1 ] with the mapping x = x i 1 2 + ξ Δ x i , y = y j 1 2 + η Δ y j , with ( ξ , η ) [ 0 , 1 ] . The way to perform, described in detail in several references, see, for instance [24], is briefly explained here.
First stage: reconstruction along x-direction. For each control volume I i j , the expressions of the reconstruction polynomials corresponding to each candidate stencil are
w h s , x ( x , t n ) = p = 0 M ψ p ( ξ ) w ^ i j , p n , s , S i j s , x ,
where ψ p are the basis interpolation functions, usually Lagrange or Legendre ones, while w ^ i j , p n , s are the coefficients of each polynomial. Integral conservation on all of the control volumes conforming each stencil yields
1 Δ x e x e 1 / 2 x e + 1 / 2 p = 0 M ψ p ( ξ ( x ) ) w ^ i j , p n , s d x = u ¯ e j n , I e j S i j s , x .
Now, we can build a data-dependent nonlinear combination of the coefficients of the polynomials for each particular stencil to yield
w h x ( x , t n ) = p = 0 M ψ p ( ξ ) w ^ i j , p n , w i t h w ^ i j , p n = s = 1 N s ω s w ^ i j , p n , s ,
where we use the nonlinear weights ω s = ω s ˜ k ω ˜ k , with ω s ˜ = λ s ( σ s + ϵ ) r , where ϵ is introduced in order to avoid division by zero. We take here ϵ = 10 20 and r = 3 , for instance, although other values can be used. The oscillation indicators are given by
σ s = p = 1 M m = 1 M Ω p m w ^ i j , p n , s w ^ i j , m n , s ,
which require the computation of the oscillation matrix
Ω p m = α = 1 M 0 1 α ψ p ( ξ ) ξ α · α ψ m ( ξ ) ξ α d ξ .
Second stage: reconstruction along y-direction. In order to obtain the reconstruction polynomial in y-direction we repeat the process followed in x-direction for each particular degree of freedom w ^ i j , p n . Subsequently, we have the expression
w h s , y ( x , t n ) = q = 0 M p = 0 M ψ p ( ξ ) ψ q ( η ) w ^ i j , p q n , s .
We now apply integral conservation in y-direction for each particular degree of freedom in x-direction, for all the control volumes of the stencil in y-direction, which is S i j s , y , to yield
1 Δ y e y e 1 / 2 y e + 1 / 2 q = 0 M ψ q ( η ( y ) ) w ^ i j , p q n , s d y = w ^ i e , p n , I i e S i j s , y ,
and perform the nonlinear combination
w h y ( x , y , t n ) = q = 0 M p = 0 M ψ p ( ξ ) ψ q ( η ) w ^ i j , p q n w i t h w ^ i j , p q n = s = 1 N s ω s w ^ i j , p q n , s .
The final expression of the reconstruction polynomial will read
w i j ( ξ , η , t n ) = k = 1 M + 1 l = 1 M + 1 w ^ i j k , l ( t n ) ψ k ( ξ ) ψ l ( η ) .
Concerning the source term, we approximate the integrals appearing in (23) while using appropriate Gaussian quadrature formulas.
In the case of the WENO5 approach, three cells are used for each stencil ( r = 3 , M = 2 ) and the developed scheme is fifth order accurate in space; whilst in the case of WENO7 we use four cells for each stencil, third degree polynomials ( r = 4 , M = 3 ) and the scheme developed is seventh order accurate in space.

4.2. Time Integration

As previously stated, time integration is achieved in this work by means of a third order Runge-Kutta TVD (Total Variation Diminishing) scheme. Relevant references on this topic are [34,35,48]. We have chosen this method, since it has good accuracy properties (third order accurate), verifies the TVD property, and is easy to implement. There are many different choices of high-order accurate methods that could be used, prominent methods are the ADER (Arbitrary high order DErivative Riemann problem) schemes that were introduced in the context of hyperbolic problems in [50,51] and subsequently applied to reaction-diffusion problems in [52,53,54], which allows obtaining arbitrary order accurate schemes in a single time step, but in the current work we have resorted to the RK3-TVD approach, since the aim of this work is to evolve the numerical solution until the stationary state is attained; hence, a fast, and still high order accurate, method is desirable. A general Runge-Kutta method to solve the ODE
d u i d t = L ( u i ) ,
where the L ( · ) is the discretized operator, is written as
u i ( j ) = k = 0 j 1 α j , k u i ( k ) + Δ t β j , k L ( u i ( k ) ) , j = 1 , m u i ( 0 ) = u i n , u i n + 1 = u i ( m ) ,
being m the number of stages in the Runge-Kutta method. The coefficients in (38) are clearly non-negative α j , k 0 , β j , k 0 . As stated in [34,35], the Runge-Kutta method (38) is TVD under the CFL (Courant–Friedrichs–Lewy) condition
c C F L = min j , k α j , k β j , k ,
where c C F L is the C F L number, provided that
α j , k 0 , β j , k 0 , k = 0 j 1 α j , k = 1 .
In the optimal case c C F L = 1 . Thus, the optimal third order Runge-Kutta TVD scheme ( m = 3 ) is given by the following expressions (see [34,35]):
  • For the solution in the surface (EBM)
    u i k , 1 = u i n + Δ t l ( t n , u i n ) , u i k , 2 = 3 4 u i n + 1 4 u i k , 1 + 1 4 Δ t l ( t n + Δ t , u i k , 1 ) , u i n + 1 = 1 3 u i n + 2 3 u i k , 2 + 2 3 Δ t l ( t n + Δ t 2 , u i k , 2 ) ,
    where l ( · ) is the discrete operator, whose expression is given in (19) and u i n refers to the cell average of the solution for the 1D control volume i at time t n .
  • For the deep ocean, we have
    U i , j k , 1 = U i , j n + Δ t L ( t n , U i , j n ) , U i , j k , 2 = 3 4 U i , j n + 1 4 U i , j k , 1 + 1 4 Δ t L t n + Δ t , U i , j k , 1 , U i , j n + 1 = 1 3 U i , j n + 2 3 U i , j k , 2 + 2 3 Δ t L t n + Δ t 2 , U i , j k , 2 ,
    where L ( · ) is the discrete operator, whose expression is given in (23) and U i , j n refers to the cell average of the solution for the control 2D volume ( i , j ) at time t n .

5. Results

The aim of this section is to numerically obtain the steady state solutions of the coupled model (8) by evolving in time the numerical scheme, according to the following Algorithm 1
Algorithm 1 Marching in time up to stationary solution
1: while Q 500 do
2:  Prescribed u i 0 , t o l
3:  for i 1 to N do
4:   while d i f > t o l do
5:    Compute u i k + 1
6:     d i f f n o r m ( u i k + 1 u i k )
7:     u i k u i k + 1
8:     k k + 1
9:   end while
10:  end for
11:  Increment the value of Q
12: end while
As a consequence, the { ( Q , u ) } diagram associated will be accomplished. With regard to the value of the tolerance, t o l , taken in the numerical simulation we have used t o l = 10 8 and computing the L 2 norm: | | u i k + 1 u i k | | 2 at each iteration. The values of the physical parameters considered in this work are given in Table 1.
Concerning the advection velocity, ω , the following function has been used
ω ( x , z ) = W ( x ) = 10 ( x + 0.75 ) ( x 0.75 ) ( 0.1 + 10 x + 0.75 ) ( 0.1 + 10 x 0.75 ) ,
which represents water sinking in the vicinity of the Poles of the Earth due to a higher density when compared to the surrounding water. Because this method is explicit stability restrictions must be taken into account. The size of the time step taken is obtained according to the expression
Δ t = min α D Δ z 2 K V 1 , α D Δ x 2 1 x 2 3 / 2 K H 0 d u d x 1 , ( α D = 0.3 ) .
We now focus on getting the stationary solutions both for the EBM model in the absence of ocean and for the coupled model (EBM+ocean). Before accomplishing this task it is determined the number of stationary solutions that can be obtained depending on the value taken by the solar constant Q. Therefore, we must consider the conditions previously stated. We choose the values: S 0 = 1 ; S 1 = 5 4 ; m = 2 and M = 0.69 , A = 190.0 , and B = 2.0 , which verify the hypotheses (11) of the theorem, as can be straightforward verified
M m S 1 S 0 1.725 1.25 1 5 x 2 4 5 4 , x M 0.4 β ( u ) 0.69 10 B + A = 170 > 0 w i t h A > 0 , B > 0 .
Therefore, we obtain the following results
Q 1 = 10 B + A M S 1 ρ c = 197.10 ; Q 2 = 10 B + A M S 0 ρ c = 246.37 , Q 3 = 10 B + A m S 1 ρ c = 340.00 ; Q 4 = 10 B + A m S 0 ρ c = 425.00 .
Hence, if Q ( 0 , Q 1 ) , there is a unique stationary solution; if Q ( Q 2 , Q 3 ) there are at least three stationary solutions and if Q ( Q 4 , + ) there is a unique stationary solution. In the following subsections, these intervals are numerically verified and the stationary solutions for different values of the solar constant Q are achieved for different initial conditions. As explained before, the numerical scheme applied is based on FVM with RK3-TVD for time integration and WENO7 spatial reconstruction.

5.1. Non-Coupled Model

We are dealing with the equation of the EBM (7) neglecting the coupling term K V U n . The model is solved for 1000 different values of the solar constant Q in the interval [ 0 , 500 ] while using eight different initial conditions.
Figure 1 depicts the { Q , u } plot attained using the following initial conditions:
  • Case 1: u ( x , 0 ) = 10.0004
  • Case 2: u ( x , 0 ) = 9.9996
  • Cases 3–6: u ( x , 0 ) = 10 + c o s ( 2 j π x ) , ( j = 0 , , 3 )
  • Case 7: u ( x , 0 ) = 10 20 c o s ( 8 π x )
  • Case 8: u ( x , 0 ) = 10 s i n ( 8 π x )
The figure shows, for different values of the solar constant Q, the L 2 norm of the difference between the vectors u and u * which represents, respectively, the stationary solution (that is, temperature) and the constant stationary solution for the value Q = 0 .
The evolution problem is solved while taking 1000 different values of the solar constant Q and the initial conditions given in Cases 1–8 until the stationary solutions are accomplished. The results are displayed in Figure 1, where those intervals of Q for which there is a unique stationary solution or a multiplicity of them are clearly identified. These regions are coherent with those theoretically predicted.
In Figure 2, Figure 3 and Figure 4, several stationary solutions for three different values of the solar constant, namely Q = 100 , 250 , 400 , are displayed. We note that both Q = 100 and Q = 400 belong to regions where the stationary solution is unique. When the value of the solar constant is Q = 250 , a multiplicity of stationary solutions appears. It can be noted that eight different solutions have been obtained, although more could be achieved while using other initial states.
In addition, it can be checked, in Figure 2, that the temperature is below zero for all of the latitudes (since the x coordinate represents the sine of the latitude). The reason behind is that, for those small values of the solar constant, the Earth receives little radiation and, as a consequence, a complete glaciation process takes place. This could have happened during the glacial periods of the Earth, which took place between 110,000 and about 15,000 years ago.
On the other hand, as displayed in Figure 4, when the solar constant take values Q > Q 4 = 425 , the unique stationary solution is positive for all latitudes and, hence, there is absence of ice throughout the Earth’s surface. This is due to the high values of the incoming solar radiation.
As expressed before, the current value of the solar constant is Q 340 , which corresponds to a multiplicity of equilibrium states.
In Figure 5, it is displayed, in a x t plot, the evolution of the solutions of the EBM model towards the stationary solution for a value of the solar constant Q = 197 . It can be verified that, for this particular value of the solar constant, the stationary solution is negative for all latitudes.
In Figure 6, it is shown the evolution towards the stationary solution taking Q = 300 , which corresponds to the region of multiple stationary solutions, taking two constant different initial conditions u ( x , 0 ) < 10 °C and u ( x , 0 ) > 10 °C. We observe that, when the initial condition is a constant below 10 °C, the stationary temperatures attained are negative for all latitudes, while when taking an initial condition constant over 10 °C, the equilibrium state achieved takes values (that is, temperatures), ranging from positive to negative, depending on the latitude.

5.2. Coupled Model

The stationary solutions of the coupled model, as given in (8), are achieved by the same procedure as in the non-coupled case that is, evolving in time the evolution problem until the steady state solution is achieved, taking 1000 values of the solar constant and using the same initial conditions (Cases 1 to 8), as before. Figure 7 shows the approximation to the bifurcation diagram accomplished in this case.
We zoom in the previous figure and achieve the plot shown in Figure 8 in order to clearly identify the regions with unique or multiple solutions.
Similar regions, as in the non-coupled situation, of unique and multiple stationary solutions show up. In order to visualize some of these solutions, the results for the values of the solar constant Q = 100 , 250 , 400 are depicted in Figure 9, Figure 10 and Figure 11.
In Figure 12, a comparison of the diagrams that were obtained for both the coupled and non-coupled model is displayed. We observe that the thermal oscillation between highest and lowest values of the temperature is larger in the coupled model than in the non-coupled one, which agrees with the thermostatic effect of the ocean.
The thermostatic effect is also patent in the results shown in Figure 13, where several stationary solutions for Q = 250 are displayed.

6. Conclusions and Discussion

In this work, we have obtained, by numerical simulation, the equilibrium states and a subset of a { ( Q , u ) } bifurcation diagram, where Q is the solar constant and u is the surface temperature. The climate model considered is composed of an Energy Balance Model and a deep ocean model. In order to obtain the stationary solutions, we have evolved in time the numerical solution computed by means of a finite volume method with a third order TVD Runge-Kutta approach for time integration, and a WENO method for spatial reconstruction.
We have achieved the stationary states of the EBM, without influence of the deep ocean, and also of the coupling EBM-deep ocean. The reason to consider these situations is that they are both relevant in climatic modelling. In the case of EBMs, without an ocean effect, they are frequently used in many studies of global climate, as they allow to predict the average surface temperature of the Earth due to incoming solar radiation, emission of radiation, energy absorption or greenhouse effects. The inclusion of the ocean effect has an important influence on the temperatures distribution due to the well-known thermostatic effect of the ocean, absorbing solar radiation and releasing it. Hence, we have dealt first with the non-coupled model, which is, solving just the Energy Balance Model without influence of the deep-ocean and then with the coupled model (surface-deep ocean), obtaining the approximation of the subset of the { ( Q , u ) } diagram, which has been achieved solving the problem taking 1000 values of the solar constant Q and using several different initial conditions for every value of Q. Three different regions have been obtained in the diagram, two of them where there is a unique equilibrium solution for each value of Q and the other one where there exists a multiplicity of equilibrium states for each particular value of Q. This result agrees with a theorem put forward by Diaz, Tello [39], where the authors prove that, in the multiplicity of equilibrium states region, there are at least three stationary solutions (that is, more than three). In the present work, we have obtained eight of those solutions, which are non-constant and stable. The results of both situations (EBM without and with the effect of the deep ocean) have been similar, although with different stationary solutions in the plot { ( Q , u ) } . Actually, a comparison have been performed between both situations, observing that, when the ocean effect is considered, the thermal oscillation, which is difference between higher and lower temperatures, is lower than when its effect is not taken into account. This agrees with the well-known thermostatic effect of the ocean. It has also been pointed out that the current situation in the Earth corresponds to the value Q 340 , which is within the multiplicity region which gives rise to temperatures that oscillate from negative to positive values in the different latitudes. Additionally, it has been indicated in this work that the region with a unique negative stationary solution for each value of Q may correspond to glaciation periods.

Author Contributions

Conceptualization, A.H. and L.T.; methodology, A.H. and L.T.; software, A.H. and L.T.; validation, A.H. and L.T.; formal analysis, A.H. and L.T.; investigation, A.H. and L.T.; resources, A.H. and L.T.; data curation, A.H. and L.T.; writing—original draft preparation, A.H. and L.T.; writing—review and editing, A.H. and L.T.; visualization, A.H. and L.T.; supervision, A.H. and L.T.; project administration, A.H. and L.T.; funding acquisition, A.H. and L.T. All authors have read and agreed to the published version of the manuscript.

Funding

The research of both authors is partially supported by the projects MTM2017-85449-P and MTM2017-83391-P Ministerio de Ciencia e Innovación, Spain.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
EBMEnergy Balance Model
DOMDeep Ocean Model
TVDTotal Variation Diminishing
WENOWeighted Essentially Non Oscillatory
RKRunge-Kutta
RK3-TVDThird order Runge-Kutta TVD
FVFinite Volume
FTCSForward in Time Centred in Space
ODEOrdinary Differential Equation
PDEPartial Differential Equation
ADERArbitrary high order DErivative Riemann problem
CFLCourant-Friedrichs-Lewy

References

  1. Budyko, M.I. The effect of solar radiation variations on the climate of the Earth. Tellus 1969, 21, 611–619. [Google Scholar] [CrossRef] [Green Version]
  2. Sellers, W.D. A Global Climatic Model Based on the Energy Balance of the Earth-Atmosphere System. J. Appl. Meteorol. 1969, 8, 392–400. [Google Scholar] [CrossRef]
  3. Díaz, J.I. On the mathematical treatment of energy balance climate models. In The Mathematics of Models for Climatology and Environment; Díaz, J.I., Ed.; Springer: Berlin/Heidelberg, Germany, 1997; pp. 217–251. [Google Scholar]
  4. North, G.R. Multiple solutions in energy balance climate models. Palaeogeogr. Palaeoclimatol. Palaeoecol. 1990, 82, 225–235. [Google Scholar] [CrossRef]
  5. Hetzer, G. The structure of the principal component for semilinear diffusion equations from energy balance climate models. Houst. J. Math. 1990, 16, 203–216. [Google Scholar]
  6. Díaz, J.; Hetzer, G.; Tello, L. An energy balance climate model with hysteresis. Nonlinear Anal. Theory Methods Appl. 2006, 64, 2053–2074. [Google Scholar] [CrossRef]
  7. Ghil, M.; Childress, S. Topics in geophysical fluid dynamics atmospheric dynamics, dynamo theory, and climate dynamics. In Applied Mathematical Sciences; Springer: New York, NY, USA, 1987. [Google Scholar] [CrossRef]
  8. Watts, R.G.; Morantine, M. Rapid climatic change and the deep ocean. Clim. Chang. 1990, 16, 83–97. [Google Scholar] [CrossRef]
  9. Díaz, J.I.; Tello, L. A nonlinear parabolic problem on a Riemannian manifold without boundary arising in climatology. Collect. Math. 1999, 50, 19–51. [Google Scholar]
  10. Diaz, J.I.; Tello, L. A 2D climate energy balance model coupled with a 3D deep ocean model. Electron. J. Differ. Equ. 2007, 16, 129–135. [Google Scholar]
  11. Bermejo, R.; Carpio, J.; Diaz, J.; Tello, L. Mathematical and numerical analysis of a nonlinear diffusive climate energy balance model. Math. Comput. Model. 2009, 49, 1180–1210. [Google Scholar] [CrossRef]
  12. Bermejo, R.; Carpio, J.; Díaz, J.I.; Galán del Sastre, P. A Finite Element Algorithm of a Nonlinear Diffusive Climate Energy Balance Model. Pure Appl. Geophys. 2008, 165, 1025–1047. [Google Scholar] [CrossRef]
  13. Díaz, J.I.; Hidalgo, A.; Tello, L. Multiple solutions and numerical analysis to the dynamic and stationary models coupling a delayed energy balance model involving latent heat and discontinuous albedo with a deep ocean. Proc. R. Soc. A 2014. [Google Scholar] [CrossRef] [Green Version]
  14. Hidalgo, A.; Tello, L. A Finite Volume Scheme for Simulating the Coupling between Deep Ocean and an Atmospheric Energy Balance Model. In Modern Mathematical Tools and Techniques in Capturing Complexity; Springer: Berlin/Heidelberg, Germany, 2011; pp. 239–255. [Google Scholar] [CrossRef]
  15. LeVeque, R.J. Finite-Volume Methods for Hyperbolic Problems; Cambridge University Press: Cambridge, UK, 2002. [Google Scholar]
  16. Godunov, S. A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations. Math. Sb. 1959, 47, 271–306. [Google Scholar]
  17. Toro, E. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  18. Vázquez-Cendón, M. Solving Hyperbolic Equations with Finite Volume Methods; Springer International Publishing: Cham, Switzerland, 2015. [Google Scholar]
  19. Harten, A.; Engquist, B.; Osher, S.; Chakravarthy, S.R. Uniformly High Order Accurate Essentially Non-oscillatory Schemes, III. J. Comput. Phys. 1997, 131, 3–47. [Google Scholar] [CrossRef]
  20. Liu, X.D.; Osher, S.; Chan, T. Weighted Essentially Non-oscillatory Schemes. J. Comput. Phys. 1994, 115, 200–212. [Google Scholar] [CrossRef] [Green Version]
  21. Shu, C. Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws. In Advanced Numerical Approximation of Nonlinear Hyperbolic Equations; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  22. Balsara, D.S.; Shu, C.W. Monotonicity Preserving Weighted Essentially Non-oscillatory Schemes with Increasingly High Order of Accuracy. J. Comput. Phys. 2000, 160, 405–452. [Google Scholar] [CrossRef] [Green Version]
  23. Titarev, V.; Toro, E. Finite-volume WENO schemes for three-dimensional conservation laws. J. Comput. Phys. 2004, 201, 238–260. [Google Scholar] [CrossRef]
  24. Dumbser, M.; Zanotti, O.; Hidalgo, A.; Balsara, D.S. ADER-WENO finite volume schemes with space–time adaptive mesh refinement. J. Comput. Phys. 2013, 248, 257–286. [Google Scholar] [CrossRef] [Green Version]
  25. Dumbser, M.; Hidalgo, A.; Zanotti, O. High order space-time adaptive ADER-WENO finite volume schemes for non-conservative hyperbolic systems. Comput. Methods Appl. Mech. Eng. 2014, 268, 359–387. [Google Scholar] [CrossRef] [Green Version]
  26. Levy, D.; Puppo, G.; Russo, G. Central WENO schemes for hyperbolic systems of conservation laws. ESAIM Modél. Math. Anal. Numér. 1999, 33, 547–571. [Google Scholar] [CrossRef]
  27. Levy, D.; Nayak, S.; Shu, C.; Zhang, Y. Central WENO schemes for Hamilton-Jacobi equations on triangular meshes. SIAM J. Sci. Comput. 2006, 28, 2229–2247. [Google Scholar] [CrossRef] [Green Version]
  28. Dumbser, M.; Boscheri, W.; Semplice, M.; Russo, G. Central Weighted ENO Schemes for Hyperbolic Conservation Laws on Fixed and Moving Unstructured Meshes. SIAM J. Sci. Comput. 2017, 39, A2564–A2591. [Google Scholar] [CrossRef]
  29. Castro, M.J.; Semplice, M. Third- and fourth-order well-balanced schemes for the shallow water equations based on the CWENO reconstruction. Math. Comput. Model. 2019, 89, 304–325. [Google Scholar] [CrossRef]
  30. Baeza, A.; Bürger, R.; Mulet, P.; Zorío, D. Central WENO Schemes Through a Global Average Weight. J. Sci. Comput. 2019, 78, 499–530. [Google Scholar] [CrossRef]
  31. Baeza, A.; Bürger, R.; Mulet, P.; Zorío, D. On the Efficient Computation of Smoothness Indicators for a Class of WENO Reconstructions. J. Sci. Comput. 2019, 80, 1240–1530. [Google Scholar] [CrossRef]
  32. Baeza, A.; Bürger, R.; Mulet, P.; Zorío, D. An Efficient Third-Order WENO Scheme with Unconditionally Optimal Accuracy. SIAM J. Sci. Comput. 2020, 42, A1028–A1051. [Google Scholar] [CrossRef]
  33. Balsara, D.S.; Garain, S.; Florinski, V.; Boscheri, W. An efficient class of WENO schemes with adaptive order for unstructured meshes. J. Comput. Phys. 2020, 404, 109062. [Google Scholar] [CrossRef]
  34. Shu, C.W.; Osher, S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J. Comput. Phys. 1988, 77, 439–471. [Google Scholar] [CrossRef] [Green Version]
  35. Gottlieb, S.; Shu, C.W. Total Variation Diminishing Runge-Kutta schemes. Math. Comput. 1998, 67, 73–85. [Google Scholar] [CrossRef] [Green Version]
  36. Harten, A. High resolution schemes for hyperbolic conservation laws. J. Comput. Phys. 1983, 49, 357–393. [Google Scholar] [CrossRef] [Green Version]
  37. Cannarsa, P.; Malfitana, M.; Martínez, P. Parameter Determination for Energy Balance Models with Memory. In Mathematical Approach to Climate Change and Its Impacts; Cannarsa, P., Mansutti, D., Provenzale, A., Eds.; Springer INdAM Series; Springer: Berlin/Heidelberg, Germany, 2020. [Google Scholar] [CrossRef] [Green Version]
  38. Floridia, G. Approximate controllability for nonlinear degenerate parabolic problems with bilinear control. J. Differ. Equ. 2014, 257, 3382–3422. [Google Scholar] [CrossRef]
  39. Díaz, J.I.; Tello, L. On a climate model with a dynamic nonlinear diffusive boundary condition. Discret. Contin. Dyn. Syst. S 2008, 1, 253. [Google Scholar] [CrossRef]
  40. Díaz, J.I.; Hernández, J.; Tello, L. On the Multiplicity of Equilibrium Solutions to a Nonlinear Diffusion Equation on a Manifold Arising in Climatology. J. Math. Anal. Appl. 1997, 216, 593–613. [Google Scholar] [CrossRef] [Green Version]
  41. Hetzer, G. The shift-semiflow of a multi-valued evolution equation from climate modelling. Nonlinear Anal. Theory Methods Appl. 2001, 47, 2905–2916. [Google Scholar] [CrossRef]
  42. Hetzer, G. The number of stationary solutions for a one-dimensional Budyko-type climate model. Nonlinear Anal. Real World Appl. 2001, 2, 259–272. [Google Scholar] [CrossRef]
  43. Díaz, J.I.; Tello, L. Infinitely many stationary solutions for a simple climate model via a shooting method. Math. Methods Appl. Sci. 2002, 25, 327–334. [Google Scholar] [CrossRef]
  44. Arcoya, D.; Diaz, J.; Tello, L. S-Shaped Bifurcation Branch in a Quasilinear Multivalued Model Arising in Climatology. J. Differ. Equ. 1998, 150, 215–225. [Google Scholar] [CrossRef] [Green Version]
  45. Stone, P.H. A simplified radiative-dynamical model for the static stability of rotating atmospheres. J. Atmos. Sci. 1972, 29, 405–418. [Google Scholar] [CrossRef] [Green Version]
  46. Díaz, J.I. Mathematical analysis of some diffusive energy balance models in Climatology. In Mathematics, Climate and Environment; Díaz, J.I., Lions, J.L., Eds.; Research Notes in Applied Mathematics n° 27; Masson: Paris, France, 1993; pp. 28–56. [Google Scholar]
  47. Hidalgo, A.; Tello, L. On a climatological energy balance model with continents distribution. Discret. Continu. Dyn. Syst. A 2015, 35, 1503. [Google Scholar] [CrossRef]
  48. Jiang, G.S.; Shu, C.W. Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef] [Green Version]
  49. Cravero, I.; Puppo, G.; Semplice, M.; Visconti, G. CWENO: Uniformly accurate reconstructions for balance laws. Math. Comp. 2018, 87, 1689–1719. [Google Scholar] [CrossRef] [Green Version]
  50. Toro, E.F.; Millington, R.C.; Nejad, L.A.M. Towards Very High Order Godunov Schemes. In Godunov Methods; Toro, E.F., Ed.; Springer: Boston, MA, USA, 2001; pp. 907–940. [Google Scholar]
  51. Titarev, V.; Toro, E. ADER: Arbitrary High Order Godunov Approach. J. Sci. Comput. 2002, 17, 609–618. [Google Scholar] [CrossRef]
  52. Toro, E.F.; Hidalgo, A. ADER finite volume schemes for nonlinear reaction diffusion equations. Appl. Numer. Math. 2009, 59, 73–100. [Google Scholar] [CrossRef]
  53. Gassner, G.; Loercher, F.; Munz, C.D. A contribution to the construction of diffusion fluxes for finite volume and Discontinuous Galerkin schemes. J. Comput. Phys. 2007, 224, 1049–1063. [Google Scholar] [CrossRef]
  54. Hidalgo, A.; Tello, L.; Toro, E.F. Numerical and analytical study of an atherosclerosis inflammatory disease model. J. Math. Biol. 2014, 68, 1785–1814. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Approximated { ( Q , | | u u * | | 2 ) } diagram for the non-coupled model. u * is the constant stationary solution for the value Q = 0 .
Figure 1. Approximated { ( Q , | | u u * | | 2 ) } diagram for the non-coupled model. u * is the constant stationary solution for the value Q = 0 .
Mathematics 08 01542 g001
Figure 2. Stationary solution of the Energy Balance Model (EBM) model, without coupling with the deep ocean, for Q = 100 .
Figure 2. Stationary solution of the Energy Balance Model (EBM) model, without coupling with the deep ocean, for Q = 100 .
Mathematics 08 01542 g002
Figure 3. Several stationary solutions of the EBM model, without coupling with the deep ocean, for Q = 250 . This is the region of multiple stationary solutions.
Figure 3. Several stationary solutions of the EBM model, without coupling with the deep ocean, for Q = 250 . This is the region of multiple stationary solutions.
Mathematics 08 01542 g003
Figure 4. Stationary solution of the EBM model, without coupling with the deep ocean, for Q = 400 .
Figure 4. Stationary solution of the EBM model, without coupling with the deep ocean, for Q = 400 .
Mathematics 08 01542 g004
Figure 5. Evolution towards the stationary solution taking as solar constant Q = 197 . The stationary temperatures are negative for all latitudes.
Figure 5. Evolution towards the stationary solution taking as solar constant Q = 197 . The stationary temperatures are negative for all latitudes.
Mathematics 08 01542 g005
Figure 6. Evolution towards the stationary solution taking as solar constant Q = 300 taking an initial condition u ( x , 0 ) < 10 °C (a) and u ( x , 0 ) > 10 °C (b).
Figure 6. Evolution towards the stationary solution taking as solar constant Q = 300 taking an initial condition u ( x , 0 ) < 10 °C (a) and u ( x , 0 ) > 10 °C (b).
Mathematics 08 01542 g006
Figure 7. Approximated { ( Q , | | u u * | | 2 ) } diagram for the coupled model.
Figure 7. Approximated { ( Q , | | u u * | | 2 ) } diagram for the coupled model.
Mathematics 08 01542 g007
Figure 8. Detail of the Approximated { ( Q , | | u u * | | 2 ) } diagram for the coupled model.
Figure 8. Detail of the Approximated { ( Q , | | u u * | | 2 ) } diagram for the coupled model.
Mathematics 08 01542 g008
Figure 9. Unique stationary solution. Coupled model. Q = 100 .
Figure 9. Unique stationary solution. Coupled model. Q = 100 .
Mathematics 08 01542 g009
Figure 10. Multiple stationary solutions. Coupled model. Q = 250 .
Figure 10. Multiple stationary solutions. Coupled model. Q = 250 .
Mathematics 08 01542 g010
Figure 11. Unique stationary solution. Coupled model. Q = 400 .
Figure 11. Unique stationary solution. Coupled model. Q = 400 .
Mathematics 08 01542 g011
Figure 12. Approximated { ( Q , | | u u * | | 2 ) } diagrams for the coupled and the non-coupled model.
Figure 12. Approximated { ( Q , | | u u * | | 2 ) } diagrams for the coupled and the non-coupled model.
Mathematics 08 01542 g012
Figure 13. Several stationary solutions for Q = 250 obtained for the EBM (a) and Coupled model (b).
Figure 13. Several stationary solutions for Q = 250 obtained for the EBM (a) and Coupled model (b).
Mathematics 08 01542 g013
Table 1. Values of the parameters used in the numerical simulation. Based on those that are given in [8], but scaled to the rectangle [ 1 , 1 ] × [ 0 , 1 ] .
Table 1. Values of the parameters used in the numerical simulation. Based on those that are given in [8], but scaled to the rectangle [ 1 , 1 ] × [ 0 , 1 ] .
ParameterValue
K H 0.049
K H 0 0.555 × 10 3
A , B 190 , 2
c, ρ 1 , 1
Q340

Share and Cite

MDPI and ACS Style

Hidalgo, A.; Tello, L. Numerical Approach of the Equilibrium Solutions of a Global Climate Model. Mathematics 2020, 8, 1542. https://doi.org/10.3390/math8091542

AMA Style

Hidalgo A, Tello L. Numerical Approach of the Equilibrium Solutions of a Global Climate Model. Mathematics. 2020; 8(9):1542. https://doi.org/10.3390/math8091542

Chicago/Turabian Style

Hidalgo, Arturo, and Lourdes Tello. 2020. "Numerical Approach of the Equilibrium Solutions of a Global Climate Model" Mathematics 8, no. 9: 1542. https://doi.org/10.3390/math8091542

APA Style

Hidalgo, A., & Tello, L. (2020). Numerical Approach of the Equilibrium Solutions of a Global Climate Model. Mathematics, 8(9), 1542. https://doi.org/10.3390/math8091542

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop