Numerical Approach of the Equilibrium Solutions of a Global Climate Model

: 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 modiﬁed 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 ﬁnite volume scheme with Weighted Essentially Non-Oscillatory (WENO) reconstruction in space and third order Runge-Kutta scheme, which veriﬁes 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 .


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 (2r − 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 n i the averaged numerical solution for cell i and time t n , it is verified that TV(u n i ) ≤ TV(u n+1 i ), where TV(·) represents the so-called Total Variation, which is given by 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.

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 − div(k(x)|∇u| p−2 ∇u) 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 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 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/m 2 . 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/m 2 . 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) = Bu + 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.

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 = sin(latitude) and z = depth. The spatial domain considered is 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 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 and for the upper boundary (surface of the ocean), we have 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, Bu + 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 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) We note here that the graph verifies the hypothesis (H 1 ).

Stationary Model
We consider the stationary problem that is associated to (8), which reads 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: this result was proved: ρc, the problem (10) has a unique solution, ρc, the problem (10) has at least three solutions, mS 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.

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 with the flux and the source term given by Concerning the deep ocean, Equation (5), we can introduce the formulation where and are the fluxes in horizontal and vertical directions respectively. Additionally, we have the source term 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 n+1 i and use these values as Dirichlet boundary condition to be used in (5) in order to obtain the cell averages U n+1 i . We integrate the Equation (7) over the space-time control volumes where is the cell average of u(x, t) over the control volume I i and is the right intercell numerical flux, while is the integral average of the source term.
We focus now on the Deep Ocean Model (DOM), which is integrated over the control volume where is the cell average of the temperature of the DOM, whereas are the spatial integral average of intercell fluxes and is integral average of the source term.

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 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  1]. The way to perform, described in detail in several references, see, for instance [24], is briefly explained here.

The following coordinates transformation is introduced
First stage: reconstruction along x-direction. For each control volume I ij , the expressions of the reconstruction polynomials corresponding to each candidate stencil are where ψ p are the basis interpolation functions, usually Lagrange or Legendre ones, whileŵ n,s ij,p are the coefficients of each polynomial. Integral conservation on all of the control volumes conforming each stencil yields Now, we can build a data-dependent nonlinear combination of the coefficients of the polynomials for each particular stencil to yield 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 which require the computation of the oscillation matrix 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ŵ n ij,p . Subsequently, we have the expression 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 and perform the nonlinear combination The final expression of the reconstruction polynomial will read 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.

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 du i dt = L(u i ), (37) where the L(·) is the discretized operator, is written as 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 where c CFL is the CFL number, provided that In the optimal case c CFL = 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) where l(·) is the discrete operator, whose expression is given in (19) and u n i refers to the cell average of the solution for the 1D control volume i at time t n .
• For the deep ocean, we have where L(·) is the discrete operator, whose expression is given in (23) and U n i,j refers to the cell average of the solution for the control 2D volume (i, j) at time t n .

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 i0 , tol 3: for i ← 1 to N do 4: while di f > tol do 5: Compute u k+1 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, tol, taken in the numerical simulation we have used tol = 10 −8 and computing the L 2 norm: ||u k+1 i − u k i || 2 at each iteration. The values of the physical parameters considered in this work are given in Table 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].

Parameter
Value Concerning the advection velocity, ω, the following function has been used 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 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 Therefore, we obtain the following results 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.

Non-Coupled Model
We are dealing with the equation of the EBM (7) neglecting the coupling term K V  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 Figures 2-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.
x Temperature (ºC) 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.

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.     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.

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. 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: