An improved method of the Globally Resolved Energy Balance model by the Bayesian networks

The accurate simulation of climate is always critically important and also a challenge. This study introduces an improved method of the Globally Resolved Energy Balance (GREB) model by the Bayesian networks based on the concept of a coarse–fine model. The improved method constructs a coarse–fine structure that combines a dynamical model with a statistical model based on employing the GREB model as the global framework and utilizing a Bayesian network constructed on the interrelationships between internal climate variables of the GREB model to achieve local optimization. To objectively validate the performance and generalization of the improved method, the method is applied to the simulation of surface temperature and temperature of the atmosphere based on the 3.75^∘ × 3.75^∘ global data sets by the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) from 1985 to 2014. The results demonstrate that the improved model exhibits higher average accuracy and lower spatial differentiation than the original GREB model and is robust in long-term simulations. This approach addresses issues with the accuracy of the GREB model in local areas, which can be attributed to an overreliance on boundary and initial conditions, as well as a lack of fully usable observed data. Additionally, the model overcomes the challenge of poor robustness in statistical models due to ambiguous climate inclusions. Thus, the improved method provides a promising way to give a reliable and stable simulation of climate.

1 Introduction

As the global warming progresses, extreme events and meteorological disasters occur frequently (Grant, 2017). Thus, the simulation and prediction of climate have become an important topic in current scientific research for the conceptual understanding and development of hypotheses for climate change studies (Dommenget and Flöter, 2011; Huang et al., 2019). Climate models are mathematical models that describe the temporal evolution of climate, oceans, atmosphere, ice, and land-use processes across a spatial domain via systems of partial differential equations (Berrocal et al., 2012), which can be solved by supercomputer and are an important tool for simulating and predicting future climate change (Kay, 2020).

Generally, climate models mainly include two categories: a dynamic model and a statistical model. A dynamic model can well understand and express the dynamic process of climate by modeling various complex climate processes or interactions, but it still faces two major problems: (i) the simulation process overly relies on initial conditions and boundary conditions (Alley et al., 2019; Zhang et al., 2019; Ludescher et al., 2021); (ii) the climate model is too complicated, and its internal characteristics cannot be fully expressed (Fan et al., 2021; Zou et al., 2019; Feng et al., 2020). The Globally Resolved Energy Balance (GREB) model is a simple but representative dynamic model, which is based on energy balance theory (Dommenget and Flöter, 2011). Compared with other dynamic models, the GREB model is a relatively fast tool for the conceptual understanding and development of hypotheses for climate change studies, because it computes about one model year per second on a standard personal computer, which allows conducting sensitivity studies to external forcing within minutes to hours (Dommenget and Flöter, 2011; Dommenget, 2016; Stassen et al., 2019). However, in addition to the two main problems of dynamic models, the GREB model also faces the problem that the model does not respond well to anomalous climate change because the parameters of the GREB model are predetermined and the observed data can hardly be used to dynamically correct the model parameters (Dommenget and Flöter, 2011; Dommenget, 2016). How to solve these problems is an important research topic to improve the GREB model and further extend it to other dynamic models.

On the contrary, a statistical model, as another type of climate models, can make good use of historical observation data to dynamically modify the models from data (Feng et al., 2020) and solve the problem that dynamic climate models rely too much on initial and boundary conditions and underutilize full observation data. Therefore, it provides a possible way to solve those defects of the dynamical model by combining that with the statistical model (Chou, 1986, 2003). A Bayesian network is a statistical method which combines graph theory and probability (Cai et al., 2013, 2019; Jansen et al., 2003). The method uses a graph to express the structure relation of the variables related to the model and has the characteristics of structuring and quantifying the object relation through the causal relation among the parts of the probability computing system (Pearl, 1986), variable logic reasoning and predictive simulation can be realized, and it can use a large amount of historical observation data. As described, it is a possible way to improve the GREB model by the Bayesian networks.

The concept of a coarse–fine model provides a joint modeling approach of dynamical–statistical hybrid model that is different from the traditional use of statistical model to optimize the empirical parameters of the dynamical model. It starts from different coarse and fine granularity of the model (Akgul and Kambhamettu, 2003; Pal and Bhattacharya, 2010; Yibo et al., 2009), uses the dynamical model as a global framework and uses the statistical model to perform local optimization, and realizes the unified modeling of both. Based on this idea, this paper introduces a method for improving the GREB model by the Bayesian networks. The aim of method is to solve the problem of low model accuracy due to overreliance on boundary conditions and initial conditions and an inability to fully utilize historical observation data. The following section presents the improved method. Section 3 presents the study case and data sets to test the new improved model. Finally, we give a discussion and conclusion of the results.

2 Methods

The improved method is developed according to the following procedure. Firstly, climate variables representing different climate processes are chosen as nodes in the Bayesian networks constructed by the GREB model. And the structural relationships among different nodes are determined to establish an abstract model of the components and structural relationships of climate processes. Secondly, the selected climate variables are categorized into variable ranges based on their numerical values to form different classifications that are used to indicate a different climate state. Thirdly, the climate state simulation method is reconstructed based on the Bayesian networks and climate evolution process to achieve the simulation of the target variable climate state. Finally, the climate state simulation results obtained from the Bayesian networks are compared with the climate model simulation results from the original GREB model to get the local optimization grids, and the numerical results of the original GREB model simulation are optimized based on the comparison results. Based on the above considerations, the improved method is developed according to the following procedures (Fig. 1).

Figure 1 Overall framework of the improved method.

2.1 Structural relationship among climate variables

Based on the energy balance, the GREB model can simulate the main characteristics and climate mean states of global warming, including seven climate processes (solar radiation, thermal radiation, hydrological cycle, sensible heat and atmospheric temperature, atmospheric circulation, sea ice, and deep ocean) and four main climatic variables (surface temperature, temperature of the atmosphere, temperature of the subsurface of the ocean, and humidity of the surface). Each of these processes is represented with strongly simplified equations. Therefore, we can abstract the structural relationship among different climate variables from the simplified equation, i.e., which climate variables control a given climate variable and which climate variables are influenced by it. This structural relationship provides the possibility to construct Bayesian networks.

2.2 Categorization of climate state

According to the theory of climate sensitivity (Annan and Hargreaves, 2006; Dommenget, 2016), climate state, indicated by a range of numerical values, can be used to replace the specific numerical values to simulate climate change and characterize the long-term trend of climate change and extreme weather conditions. And it is better suited to capture the similarity of a given climate variable across different spatial and temporal locations compared to the numerical values of specific climate variables. Therefore, it can be used to assess the similarity between the simulated results of a model and the actual results, indicating the accuracy of the simulation. This provides a simple and practical approach to evaluating the accuracy of revealing local abrupt changes in simulation results. Moreover, by simulating state rather than a specific numerical value, it is possible to significantly reduce computational effort and simulation response time. This is consistent with the primary objective of the GREB model, which is to provide a fast tool for the conceptual understanding and development of hypotheses for climate change studies (Dommenget and Flöter, 2011).

The natural breaks classification (Jenks) method is a commonly used classification method that aims to minimize intra-class variation and maximize inter-class variation. By categorizing the numeric values of climate variables into different classifications to indicate the climate state using the natural breaks classification method, it can be considered that numeric values within the same classification have less variation, representing the fact that the results of this classification of numeric values have a similar climate state.

2.3 Climate state simulation

According to the characteristics of Bayesian networks, the climate state simulation of climate variables is realized by a climate evolution process based on Bayesian networks; i.e., the climate state of unknown climate variables is inferred from the climate state of known climate variables at the same spatial locations.

2.3.1 Bayesian networks

A Bayesian network is a probabilistic model that simulates the human reasoning process, which is a combination of graph theory and probability theory, and its network topology is a directed acyclic graph, where variables are nodes and correlations or causal relationships between variables are directed edges. The dynamic evolution of Bayesian network node probabilities is controlled by conditional probabilities, and each node covers a probability distribution table under the joint distribution of the parent nodes, indicating the strength of the relationship between the nodes (Sahin et al., 2019). When the Bayesian network is constructed, given the state of any node, the probability distribution of the states of the remaining nodes can be calculated.

In the Bayesian networks, the probability of a node can be calculated in the form of probability using prior knowledge and statistical data, namely the Bayes probability (Maher, 2010). Observed sample are defined as $G=\left\{X_{\mathrm{1}}=x_{\mathrm{1}},X_{\mathrm{2}}=x_{\mathrm{2}},\mathrm{\cdots },X_n=x_n\right\}$, where X is the event and x is the event value or state. When θ is the prior probability of event X=x, ζ is prior knowledge, and P(θ|ζ) is the probability density function, then the probability $P\left(X_{n+\mathrm{1}}=x_{n+\mathrm{1}}|\mathit{\theta },\mathit{\zeta }\right)$ of the n+1 event $X_{n+\mathrm{1}}=x_{n+\mathrm{1}}$ can be obtained from the prior probability density P(θ|ζ) and the sample G through the Bayes probability. It can be calculated by the total probability formula:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}P(X_{n+\mathrm{1}}& =x_{n+\mathrm{1}}|\mathit{\theta },\mathit{\zeta })\\ & =\int P\left(X_{n+\mathrm{1}}=x_{n+\mathrm{1}}|\mathit{\theta },G,\mathit{\zeta }\right)P\left(\mathit{\theta }|,G,\mathit{\zeta }\right)\mathrm{d}\mathit{\theta }\\ & =\int \mathit{\theta }P\left(\mathit{\theta }|,G,\mathit{\zeta }\right)\mathrm{d}\mathit{\theta }.\end{array}\end{array}$$

Based on the Bayes equation, the posterior probability $P\left(\mathit{\theta }|,G,\mathit{\zeta }\right)$ is denoted as

$$\begin{array}{}\text{(2)}& P\left(\mathit{\theta }|,G,\mathit{\zeta }\right)={\displaystyle \frac{P\left(\mathit{\theta }|\mathit{\zeta }\right)P\left(G|\mathit{\theta }\text{,}\mathit{\zeta }\right)}{P\left(G|\mathit{\zeta }\right)}},\end{array}$$

where G is the given sample, ζ is the prior probability of G, and θ is the posterior probability of G.

2.3.2 Climate evolution process based on Bayesian networks

In a climatic process composed of several climatic variables, there is an association relationship between climatic variables. These climatic variables are regarded as network nodes, and the association relations between climatic variables are taken as directed edges. The association relationship between nodes is represented by the graph model, and the action intensity of the association relationship is described quantitatively by the conditional probability table. Using the characteristics of Bayesian networks, the attribute feature state of nodes is inferred by probability to realize the expression and simulation of the attribute feature state of geographical variables.

A climate process $M_t=\left\{X\left(m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_i\right)|m_{\mathrm{1}t},m_{\mathrm{2}t},\mathrm{\dots },m_{it}\right\}$ is composed with i climate variables, $m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_i$, and $X\left(m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_i\right)$ is the structural relationship among the variables. Suppose that the climate variable m_i has j states, then the states set of m_i is $\left\{W_{m_{i\mathrm{1}}},W_{m_{i\mathrm{2}}},\mathrm{\dots },W_{m_{ij}}\right\}$. The climate process is described by a Bayesian network $B=\left(S,X\right)$, where S is a directed acyclic graph composed of nodes; X is the node set of the graph, that is, climate variables $m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_i$. Nodes are connected by directed edges to represent the relationship between climate variables. Each node has an independent conditional probability table, which represents the probability distribution under the joint distribution of its parent nodes. Assume that a climate m_i has one or more parent nodes $m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_e\left(e\le i-\mathrm{1}\right)$ and states $d_{\mathrm{1}},d_{\mathrm{2}},\mathrm{\dots }{d}_e$, and it can be denoted as $m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_e\to m_i$. Under the parent node of all possible states, the conditional probability table composed of the set of state probabilities of m_i is as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& B_{m_i}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}\\ & =\left\{(W_{m_{i\mathrm{1}}},P_{W_{m_{i\mathrm{1}}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}),\right.\\ & (W_{m_{i\mathrm{2}}},P_{W_{m_{i\mathrm{2}}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}),\mathrm{\dots },\\ & \left.(W_{m_{ij}},P_{W_{m_{ij}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}})\right\}\\ & \left(r\mathrm{1}=\mathrm{1},\mathrm{2},\mathrm{\dots },d_{\mathrm{1}}\right),\left(r\mathrm{2}=\mathrm{1},\mathrm{2},\mathrm{\dots },d_{\mathrm{2}}\right),\mathrm{\dots },\\ & \left(re=\mathrm{1},\mathrm{2},\mathrm{\dots },d_e\right),\end{array}\end{array}$$

where $B_{m_i}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}$ is a conditional probability table of climate variables m_i; $W_{m_{ij}}$ is the jth characteristic state of climate variables m_i; and $P_{W_{m_{ij}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}$ is the probability of climate variable m_i that corresponds to the jth state under the $r\mathrm{1},r\mathrm{2},\mathrm{\dots }re$ characteristic state corresponding to the parent node $m_{\mathrm{1}},m_{\mathrm{2}},\mathrm{\dots },m_e$ expression set. The probability set of climate variable m_i at t moment can be denoted as $C_{m_{it}}$:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& C_{m_{it}}=\\ & \left\{(W_{m_{i\mathrm{1}}},P_{W_{m_{i\text{1}}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}),\right.\\ & (W_{m_{\text{i2}}},P_{W_{m_{i\mathrm{2}}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}}),\mathrm{\dots },\\ & \left.(W_{m_{ij}},P_{W_{m_{ij}}}^{W_{m_{\mathrm{1}r\mathrm{1}}},W_{m_{\mathrm{2}r\mathrm{2}}},\mathrm{\dots },W_{m_{ere}}})\right\}.\end{array}\end{array}$$

The conditional probability table of each node can be calculated by Eq. (2) using training data.

2.4 Local optimization

The numerical results simulated by the original GREB model are compared with the climate state results simulated by the Bayesian networks, and the grids where the numerical result simulated by the original GREB model are not in the range of the climate state simulated by the Bayesian networks are used as grids to be optimized.

According to the Third Law of Geography (Zhu et al., 2018), the more similar the geographic environment, the more similar the geographic target characteristics are. Therefore, for an unknown climate variable at a certain spatial and temporal location, the numerical values of other known climate variables at that spatial and temporal location can be used to infer. Accordingly, we propose that for an unknown climate variable, the position of its specific value in the range of its classification is related to the position of the specific value of the known climate variable in the range of its classification at the same spatial and temporal location. For a climate variable containing n relevant control variables, the numerical results are calculated as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{E}_{\mathrm{value}}^x& =S_{\mathrm{lower}\mathrm{limit}}^E+{\displaystyle \frac{\mathrm{1}}{n}}\left(S_{\mathrm{upper}\mathrm{limit}}^E-S_{\mathrm{lower}\mathrm{limit}}^E\right)\\ & \cdot \sum {\displaystyle \frac{E_{\mathrm{value}}^i-S_{\mathrm{lower}\mathrm{limit}}^i}{S_{\mathrm{upper}\mathrm{limit}}^i-S_{\mathrm{lower}\mathrm{limit}}^i}},\end{array}\end{array}$$

where $E_{\mathrm{value}}^x$ represents an unknown climate variable, $S_{\text{lower limit}}^E$ represents the lower limit of the range of classification in which the unknown climate variables are simulated by Bayesian networks, $S_{\text{upper limit}}^E$ represents the lower limit of the range of classification in which the unknown climate variables are simulated by Bayesian networks, n represents the number of known climate variables associated with the unknown variables in the Bayesian networks, $E_{\mathrm{value}}^i$ represents the actual value of the ith known climate variable, $S_{\text{lower limit}}^i$ represents the lower limit of the range of classification in which the ith known climate variables are in the simulation process by Bayesian networks, and $S_{\text{upper limit}}^i$ represents the upper limit of the range of classification in which the ith known climate variables are in the simulation process by Bayesian networks.

According to the above method, we can improve the accuracy of the model by comparing the climate state, identifying the grid to be optimized, and recalculating the values simulated by the original GREB model within the grid. In this way, the improved model with coarse–fine structure constructed with the GREB model as the global framework and the Bayesian networks as the local optimization can better reflect the localized abrupt changes in the climate process and achieve the purpose of improving the GREB model.

3 Case study

In order to demonstrate the accuracy of the improved model in simulating climate variables and to verify its reliability, surface temperatures and temperature of the atmosphere from the GREB model were selected for simulation objects. The simulation of these two climate variables includes most of the climate processes of the GREB model and can reflect the complex coupling process and climate change characteristics of the GREB model.

3.1 Data description

In this paper, data produced by the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) are used as the experimental data to evaluate the improved model. The data sets include surface temperature (T_surf), temperature of the atmosphere (T_atmos), solar radiation (F_solar), total cloud cover (CLD), water vapor (q_air), temperature of the subsurface ocean (T_ocean), and wind speed (u) stored as a 3.75^∘ × 3.75^∘ (latitude × longitude) grid NC data from 1985 to 2014. In order to facilitate calculation and comparative analysis, all climate data are preprocessed. Firstly, the downloaded climate data are removed from the outliers so that the data are calculated to avoid too large or too small results; secondly, the grid data are resampled, and the resampling method is bilinear interpolation. The bilinear interpolation method is used to interpolate the climate data, which not only fills the null values, but also unifies the scale size of the data. Finally, considering that changes in climate variables are usually seasonally related, climate data from 1985 to 2014 were processed as quarterly averages, where January, February, and March comprised the first quarter; April, May, and June formed the second quarter; July, August, and September constituted a third quarter; and October, November, and December comprised the fourth quarter.

3.2 Structural relationship among climate variables and climate state

The process of simulating the surface temperature includes solar radiation, thermal radiation, sensible heat and atmospheric temperature, and deep ocean (Dommenget and Flöter, 2011). The main heat source of the surface temperature is solar radiation: some is absorbed by the surface temperature; the other part is reflected by the surface temperature, and part of the heat on the surface temperature is transferred in the atmosphere; and some of it is transferred to the ocean below the surface. Each climate variable in this scene can be expressed by a highly simplified equation, which follows the surface temperature tendency equation as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\gamma }}_{\mathrm{surf}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{surf}}}{\mathrm{d}t}}=F_{\mathrm{solor}}+F{a}_{\mathrm{thermal}}+F_{\mathrm{latent}}+F_{\mathrm{sense}}+F_{\mathrm{ocean}},\end{array}$$

where T_surf is surface temperature, γ_surf is surface heat capacity, F_solar is the incoming solar radiation, F_thermal is the net thermal radiation, F_latent is the cooling by latent heat from surface evaporation of water, F_sense is the turbulent heat exchange with the atmosphere, and F_ocean is the heat exchange with the deeper subsurface ocean. The subprocesses of surface temperature are modeled as follows:

$$\begin{array}{}\text{(7)}& \left\{\begin{array}{l}{F}_{\mathrm{solar}}=\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{clouds}}\right)\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{surf}}\right)S_{\mathrm{0}}\cdot r\left(\mathit{\varphi },t_{\mathrm{julian}}\right)\\ F_{\mathrm{thermal}}=-\mathit{\sigma }{T}_{\mathrm{surf}}^{\mathrm{4}}+{\mathit{\epsilon }}_{\mathrm{atmos}}\mathit{\sigma }{T}_{\mathrm{atmos}-\mathrm{rad}}^{\mathrm{4}}\\ {\mathit{\epsilon }}_{\mathrm{atmos}}=\frac{p{e}_{\mathrm{8}}-\text{CLD}}{p{e}_{\mathrm{9}}}\cdot \left({\mathit{\epsilon }}_{\mathrm{0}}-p{e}_{\mathrm{10}}\right)+p{e}_{\mathrm{10}}\\ F_{\mathrm{latent}}=L\cdot {\mathit{\rho }}_{\mathrm{air}}\cdot C_w\cdot \left|{\mathit{u}}_{*}\right|\cdot {\mathit{\upsilon }}_{\mathrm{soil}}\cdot \left(q_{\mathrm{air}}-q_{\mathrm{sat}}\right)\\ F_{\mathrm{sense}}=c_{\mathrm{atmos}}\cdot \left(T_{\mathrm{atmos}}-T_{\mathrm{surf}}\right)\\ F_{\mathrm{ocean}}=F{o}_{\mathrm{sense}}+{\mathit{\gamma }}_{\mathrm{surf}}\cdot \mathrm{\Delta }{T}_{\mathrm{entrain}}\\ F{o}_{\mathrm{sense}}=c_{\mathrm{ocean}}\cdot \left(T_{\mathrm{ocean}}-T_{\mathrm{surf}}\right)\end{array}\right.,\end{array}$$

where F_solar is the incoming solar radiation, α_clouds is the fraction of the incoming solar radiation reflected by clouds, α_surf is the fraction of the incoming solar radiation reflected by the surface, S₀ is the solar constant, α_surf is the fraction of the incoming solar radiation reflected by the surface, r is the 24 h mean fraction reaching a normal surface area on top of the atmosphere, ϕ is the function of latitude, t_julian is the Julian day of the calendar year, F_thermal is the net thermal radiation, T_surf is surface temperature, ε_atmos is the effective emissivity, T_atmos−rad is the temperature defined in the context of the atmospheric temperature, CLD is the total cloud cover, ε₀ is the emissivity without considering clouds first, pe_* is the parameters, F_latent is the cooling by latent heat from surface evaporation of water, L is the constant parameters of the latent heat of evaporation and condensation of water, ρ_air is the density of air, C_w is the transfer coefficient, |u_*| is the wind speed, υ_soil is the bulk formula extended by a surface wetness fraction, q_air is the actual surface air layer humidity, q_sat is the saturation surface air layer specific humidity, F_sense is the turbulent heat exchange with the atmosphere, c_atmos is the coupling constant, T_atmos is temperature of the atmosphere, F_ocean is the heat exchange with the deeper subsurface ocean, Fo_sense is the turbulent mixing between the two ocean layers, ΔT_entrain is the heat exchange with the surface ocean layer due to decreasing of the mixed layer depth, c_ocean is the coupling constant, and T_ocean is the temperature of the subsurface ocean.

In the process of simulating the temperature of the atmosphere by the GREB model (Dommenget and Flöter, 2011), the temperature of the atmosphere is not only related to the thermal radiation reflected from the surface, but also related to the sensible heat exchange with the surface and latent heat release by condensation of atmospheric water vapor. Each climate variable in this process can be expressed by a highly simplified equation, which follows the temperature of the atmosphere tendency equation as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{\gamma }}_{\mathrm{atmos}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{atmos}}}{\mathrm{d}t}}& +F_{\mathrm{sense}}=F{}_{\mathrm{thermal}}+Q_{\mathrm{latent}}\\ & +{\mathit{\gamma }}_{\mathrm{atmos}}\left(k\cdot {\mathrm{\nabla }}^{\mathrm{2}}{T}_{\mathrm{atmos}}-\mathit{u}\cdot \mathrm{\nabla }{T}_{\mathrm{atmos}}\right),\end{array}\end{array}$$

where T_atmos is temperature of the atmosphere, γ_atmos is atmospheric heat capacity, F_sense is the sensible heat exchange with the surface, Fa_thermal is net thermal radiation of the atmosphere, Q_latent is the latent heat release by condensation of atmospheric water vapor, and u is the wind speed. The subprocesses of temperature of the atmosphere are modeled as follows:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}{F}_{\mathrm{sense}}=c_{\mathrm{atmos}}\cdot \left(T_{\mathrm{atmos}}-T_{\mathrm{surf}}\right)\\ F{a}_{\mathrm{thermal}}={\mathit{\epsilon }}_{\mathrm{atmos}}\mathit{\sigma }\cdot T_{\mathrm{surf}}^{\mathrm{4}}-\mathrm{2}{\mathit{\epsilon }}_{\mathrm{atmos}}\mathit{\sigma }\cdot T_{\mathrm{atmos}-\mathrm{rad}}^{\mathrm{4}}\\ Q_{\mathrm{latent}}=-\mathrm{2.6736}\times {\mathrm{10}}^{\mathrm{3}}\left[\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}\right]\cdot \mathrm{\Delta }{q}_{\mathrm{precip}}\cdot L\end{array}\right.,\end{array}$$

where F_sense is the turbulent heat exchange with the atmosphere, c_atmos is the coupling constant, T_atmos is temperature of the atmosphere, Fa_thermal is net thermal radiation of the atmosphere, T_surf is surface temperature, ε_atmos is the effective emissivity, T_atmos−rad is the temperature defined in the context of the atmospheric temperature, Q_latent is the latent heat release by condensation of atmospheric water vapor, Δq_precip is the condensation or precipitation, and L is the constant parameters of the latent heat of evaporation and condensation of water.

For different climate processes, the climate subprocesses and relationship structures are different. Therefore, the selection of nodes in each climate process will also be different. Not only is the selection of appropriate variables as nodes very important, but also the number of nodes will directly affect the simulation of the final climate average state. In order to simplify the complex climate evolution process and facilitate calculation, four to six climate variables are selected as key nodes in each climate process, and the variable climate state in each process is simulated by these nodes.

Through the trend equations (Eqs. 6 and 7) in the processes of surface temperature, the relation equation of climate variables can be simplified:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{l}{T}_{\mathrm{surf}}=f(F_{\mathrm{solar}},T_{\mathrm{ocean}},q_{\mathrm{air}},\text{CLD})\\ T_{\mathrm{ocean}}=f\left(F_{\mathrm{solar}}\right)\\ q_{\mathrm{air}}=f\left(F_{\mathrm{solar}}\right)\end{array}\right.,\end{array}$$

where T_surf is the surface temperature; F_solar is solar radiation; T_ocean is the temperature of the subsurface ocean; q_air is the actual surface air layer humidity, i.e., water vapor content; and CLD is the total cloud cover. That is, the surface temperature, solar radiation, temperature of the subsurface ocean, total cloud cover, and water vapor content can be selected as the key nodes of the surface temperature process.

Through the trend equation (Eqs. 8 and 9) in the processes of the temperature of the atmosphere, the relation equation of climate variables can be simplified:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}{T}_{\mathrm{atmos}}=f(\mathit{u},q_{\mathrm{air}},\text{CLD})\\ q_{\mathrm{air}}=f\left(\mathit{u}\right)\end{array}\right.,\end{array}$$

where T_atmos is the temperature of the atmosphere, u is the wind speed, CLD is the total cloud cover, and q_air is the water vapor content. That is, the temperature of the atmosphere, wind speed, total cloud cover, and water vapor content can be selected as the key nodes of the temperature of the atmosphere process.

According to Eq. (10), in the surface temperature process, the surface temperature (T_surf) is controlled by solar radiation (F_solar), total cloud cover (CLD), water vapor (q_air), and temperature of the subsurface ocean (T_ocean). The temperature of the subsurface ocean (T_ocean) and water vapor (q_air) are controlled by solar radiation (F_solar). For the above relationship, the Bayesian network structure in the surface temperature process can be constructed (Fig. 2a). According to Eq. (11), in the temperature of the atmosphere process, the temperature of the atmosphere (T_atmos) is controlled by cloud cover (CLD), water vapor (q_air), and wind speed (u). And water vapor (q_air) is controlled by wind speed (u). For the above relationship, the Bayesian network structure in the temperature of the atmosphere process can be constructed (Fig. 2b).

Figure 2 Structural relationship among climate variables in the different simulation processes. (a) Surface temperature process; (b) temperature of the atmosphere process.

The climate state of the variables in the above climate processes was performed using the natural breaks classification method. The climate variable data are categorized into five, seven, and nine different classifications to indicate different climate states to test the improved model and verify the effect of the classification number of climate variable data on the simulation results. Detailed schemes are shown in Appendix Tables B1, B2, and B3.

3.3 Climate state simulation

Surface temperature and temperature of the atmosphere are considered as a simulation object and other climate variables as known objects, and historical data are used to calculate the conditional probability tables of each nodes through the Bayesian network structure with Eq. (4). Among them, the training data are the 10-year historical data from 1985 to 1994.

In each simulation process, there are two training methods for the simulated object. The first is to train a conditional probability table using the data in all the grids, and the second is to use the conditional probability table to simulate the states of all grids. The conditional probability table obtained by this training method can reflect the numerical characteristic relationship between climate variables in the whole region. However, it cannot show the distribution pattern of the characteristics of the simulated state in space. The second is to train the data in each grid separately. Because the state grading data in each grid is different, the conditional probability table of the simulated object trained in each grid is also different, and a total of 96×48 conditional probability tables are obtained. The conditional probability tables obtained by this training method can accurately reflect the different numerical characteristic relationship between the simulated object and the known object in different regions. However, due to the training of more conditional probability tables, the running time of this training method will be a little longer. Considering the great differences in the pattern of climate evolution in different regions, this paper uses the second data training method in state simulation, which first divides the whole world into 96×48 grids and then uses the data in each grid to train the conditional probability tables of the grid. After the training is completed, the data of the known climate variables will be used to simulate the unknown climate variables from 1995 to 2014. The simulation results are shown in Figs. 3 and 4.

Figure 3 Climate state simulation results of the surface temperature by Bayesian networks for 80 seasons (the period 1995–2014). (a) Categorized into five classifications; (b) categorized into seven classifications; (c) categorized into nine classifications.

Figure 4 Climate state simulation results of the temperature of the atmosphere by Bayesian networks for 80 seasons (the period 1995–2014). (a) Categorized into five classifications; (b) categorized into seven classifications; (c) categorized into nine classifications.

Figure 3 shows the climate state of the quarterly average of surface temperature from 1995–2014. The simulation results under different classifications all clearly show the global quarterly average surface temperature distribution with latitudinal variations. The surface temperature starts from the Equator and decreases with the increase in latitude, so the temperature in the North and South Pole is the lowest. The climate state distribution of surface temperature is basically in line with the real world. Different from the simulation result of surface temperature, the quarterly average temperature of the atmosphere rises from the Equator and increases with the increase in latitude in Fig. 4, which is also basically in line with the real world. The tropospheric height of the poles is lower, and the tropospheric height of the Equator is higher – a phenomenon which leads to the result that the temperature of the troposphere at the same height is higher in the poles.

3.4 Local optimization

After the climate state simulation of Bayesian networks, the numerical results simulated by the original GREB model are compared with the climate state results simulated by the Bayesian networks, and the grids where the numerical results simulated by the original GREB model are not in the range of the climate state simulated by the Bayesian networks are used as grids to be optimized. The GREB model uses the model code sourced from the Monash Simple Climate Model (MSCM) laboratory repository for its implementation, and we run the code in the FORTRAN language.

Based on the optimized area of surface temperature simulations and temperature of the atmosphere obtained from the climate state accuracy comparison (see Appendix A for details), the original GREB model simulation results of the grid to be optimized in the optimized area are recalculated according to Eq. (5).

In terms of surface temperature simulation, the original GREB model at low latitudes shows high state accuracy, so a local optimization scheme is used only for the middle and high latitudes. The empirical parameter for the optimization range of the surface temperature simulation has been determined to be 90 to 30^∘ N and 30 to 90^∘ S. The quarterly average surface average temperature simulated by the improved model for the period 1994–2015 is presented in Fig. 5. In terms of the temperature of the atmosphere simulation, a spatially global optimization approach has been chosen, owing to the higher global state accuracy of the Bayesian networks. The quarterly average temperature of the atmosphere simulated by the improved model for the period 1994–2015 is presented in Fig. 6. The details in Figs. 5 and 6 show that the optimization data results are well characterized by localized abrupt changes, which means that the improved model is able to effectively address the inadequate response of the original GREB model to localized abrupt changes.

Figure 5 Quarterly average surface temperature for 80 seasons (the period 1995–2014). (a) Simulated by the improved GREB model based on Bayesian networks under five classifications; (b) simulated by the improved GREB model based on Bayesian networks under seven classifications; (c) simulated by the improved GREB model based on Bayesian networks under nine classifications.

Figure 6 Quarterly average temperature of the atmosphere for 80 seasons (the period 1995–2014). (a) Simulated by the improved GREB model based on Bayesian networks under five classifications; (b) simulated by the improved GREB model based on Bayesian networks under seven classifications; (c) simulated by the improved GREB model based on Bayesian networks under nine classifications.

3.5 Evaluation of improved model

In order to evaluate the simulation accuracy of the improved model (the optimized GREB model based on Bayesian networks of climate state), the root mean square error (RMSE) between the simulated and actual values is defined to evaluate the model:

$$\begin{array}{}\text{(12)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum (S_i-A_i{)}^{\mathrm{2}}},\end{array}$$

where S_i represents the simulated value and A_i represents the actual value; when analyzed spatially n represents the length of time, whereas when analyzed temporally n represents the number of grids in space. The accuracy of original GREB model simulation result was used as a comparison.

The mean values of the RMSE between the simulated results of the original GREB model, as well as the improved models based on five, seven, and nine classifications, and the observed values for surface temperature were 13.26, 8.66, 8.85, and 9.81, respectively. For the temperature of the atmosphere, the corresponding mean values of the RMSE were 72.19, 22.77, 20.12, and 17.76, respectively. This result shows that the improved method significantly reduces the RMSE of the simulation; i.e., it improves the simulation accuracy. However, there are also significant differences between the surface temperature and the temperature of the atmosphere. There is no significant relationship between the RMSE and classification in the simulation of surface temperature, while the RMSE decreases with increasing classification in the simulation of temperature of the atmosphere.

Figure 7 depicts the spatial distribution of the RMSE between the simulated surface temperature and the observed values for the original GREB model (Fig. 7a) and the improved models based on five, seven, and nine classifications (Fig. 7b–d). The comparison shows that the improved model significantly improves the simulation accuracy of the surface temperature in the polar regions. Figure 8 depicts the spatial distribution of the RMSE between the simulated temperature of the atmosphere and the observed values for the original GREB model (Fig. 8a) and the improved models based on five, seven, and nine classifications (Fig. 8b–d). The comparison shows that the improved model significantly improves the simulation accuracy of temperature of the atmosphere at mid- and low-latitude regions. This is also well verified by the variation curve of the RMSE along the latitude direction shown in Fig. 9.

Figure 7 Spatial distribution of the RMSE for 80 seasons of the simulated surface temperature. (a) Simulated by the original GREB model; (b) simulated by the optimized GREB model based on Bayesian networks under seven classifications; (c) simulated by the optimized GREB model based on Bayesian networks under five classifications; (d) simulated by the optimized GREB model based on Bayesian networks under nine classifications.

Figure 8 Spatial distribution of the RMSE for 80 seasons of the simulated temperature of the atmosphere. (a) Simulated by the original GREB model; (b) simulated by the optimized GREB model based on Bayesian networks under five classifications; (c) simulated by the optimized GREB model based on Bayesian networks under seven classifications; (d) simulated by the optimized GREB model based on Bayesian networks under nine classifications.

Figure 9 Mean of the spatial distribution of the RMSE in the latitudinal direction; GREB represents the results of the original GREB model simulation, IMPM-5 represents the simulation results of the improved model under five classifications, IMPM-7 represents the simulation results of the improved model under seven classifications, and IMPM-9 represents the simulation results of the improved model under nine classifications. (a) Comparison results of the simulated temperature; (b) comparison results of the simulated temperature of the atmosphere.

Figure 10 shows the quarterly variability and trends of the RMSE between 1995 and 2014 for both the surface temperature and temperature of the atmosphere. The comparison results demonstrate that the improved model significantly reduces the RMSE and exhibits temporal stability, indicating the robustness of the improved model. Moreover, the RMSE curves of the improved models exhibit the same seasonal cycle as the original GREB model, with the smallest RMSE occurring in the fourth quarter and the largest in the third quarter. This seasonal pattern can be attributed to the fact that the improved model is based on the modeling of the climate variable relationship within the GREB model, thus exhibiting similar temporal variation characteristics to those of the GREB model, which reflects the coarse–fine structure of improved model with the original GREB model as the global framework. The RMSE trends over time demonstrate that the improved model is temporally stable, and its accuracy does not deviate over time. This renders the improved model suitable for simulating surface temperature and temperature of the atmosphere over long time series.

Figure 10 Quarterly variability and trends of the RMSE between 1995 and 2014; GREB represents the results of the original GREB model simulation, IMPM-5 represents the results of the improved model under five classifications, IMPM-7 represents the results of the improved model under seven classifications, IMPM-9 represents the results of the improved model under nine classifications. (a) Comparison results of the simulated surface average temperature; (b) comparison results of the simulated atmospheric average temperature.

4 Conclusions and discussions

In this study, we introduced a coarse–fine structure to improve the GREB model based on Bayesian networks. The improved model uses the GREB model as the basis of the global simulation framework and uses the Bayesian networks to do local optimization. By introducing a Bayesian networks, the results of the original GREB model are quickly evaluated with the climate state as the evaluation index, the local optimization region is confirmed, and the simulation results of the GREB model within the optimization region are recalculated, which improves the model accuracy significantly.

The improved model was evaluated by two cases: surface temperature and temperature of the atmosphere. The simulation results of the improved model show that the improved model has higher average accuracy and lower spatial variability compared to the original GREB model. This means that the improved model has better applicability and stability on a global scale. Meanwhile, on the timescale, the model maintains good robustness and does not suffer from the problem of accuracy divergence of traditional statistical models because the improved model uses the GREB model as the basic global framework. The results of the two study cases not only demonstrate that the improved method can also be used for the simulation of other climate variables within the GREB model, but they also reveal the construction of coarse–fine models through a combination of dynamical and statistical methods as a potential means of improving climate simulation and prediction. This improved approach can overcome the shortcomings of a single dynamical model that cannot accurately describe many nonlinear processes in the climate system and can be applied to other dynamical models. In terms of development, the improved methods for improving climate dynamical models by statistical methods show great possibilities for improving the accuracy of climate predictions.

In addition to the improved model with improved accuracy, the concept of evaluation through climate state introduced during the construction of the coarse–fine model can be well in studies on climate sensitivity (Dommenget, 2016; Kutzbach et al., 2013), extreme weather (Bellprat and Doblas-Reyes, 2016; Chen et al., 2018), and climate threshold (Mahlstein et al., 2015; Vogel et al., 2020). In particular, climate anomalies, as manifested by the climate state, can serve as effective indicators for tracking climate change trends. However, this improved method still has some shortcomings.

1. The scientific problem of categorization of climate variable attribute features. In this paper, the climate state of each variables is indicated by the classifications categorized by the natural breaks classification method according to the data characteristics and statistical regularities of the climate model, but this classification method changes with the data, and the data-based classification model may not be consistent with the actual climate evolution pattern. Therefore, the following studies can discuss related issues and choose the appropriate feature classification criteria to achieve a balance between different simulation.

2. Balance of accuracy and resolution. If the actual numerical values rather than states are used as the calculation parameters, a higher resolution will be obtained, and of course the training data for each case will be reduced, which leads to the loss of accuracy. How to achieve the balance of accuracy and resolution will be an important issue.

3. Applicability of climate evolution models based on Bayesian networks. Stable conditional probability tables can be trained with historical climate data to simulate climate state, but conditional probability tables cannot change over time and cannot be adapted to time-sensitive climate models. The following study can extend the applicability of the method by dynamically training Bayesian networks on climate data.

Appendix A: State accuracy comparison

In order to verify the reliability of the simulated climate state using the Bayesian networks and to provide a basis for guiding the optimization of the GREB local simulation result, the state accuracy (dimensionless) was used to evaluate the reliability of the simulated climate state, which is expressed as

$$\begin{array}{}\text{(A1)}& \text{state accuracy}={\displaystyle \frac{n}{N}},\end{array}$$

where n represents the number of time series in which the simulated state value of a grid is the same as the actual state value in the time series, which in this case refers to the number of seasons; and N represents the total number of time series. State accuracy means the same proportion of the simulated and actual states in the same grid.

The numerical results simulated by the original GREB model are also transformed into the climate state by the natural breaks classification method for comparative evaluation. The state accuracy averaged over different processes from 1985 to 1994 is shown in Fig. A1, and the state accuracy of the surface temperature and the temperature of the atmosphere are shown in Figs. A2 and A3. Overall, the comparison results (Fig. A1) show that the Bayesian networks have a higher simulation state accuracy in both the surface temperature and the temperature of the atmosphere. This higher state accuracy indicates that the Bayesian networks simulate climate state better than the GREB model, which provides a basis for evaluating the GREB model simulations with Bayesian network simulation results. Intrinsic to this result is the fact that since more observations are involved in the simulation process (in the construction of conditional probability tables) in Bayesian network simulations, this allows the Bayesian network response to localized abrupt changes in climate to be more pronounced.

When it comes to the number of classifications, the total number of data remains unchanged, and as the number of classifications increases, the number of training data per classification decreases, which results in a decrease in the accuracy of the simulations of the two methods. This implies that the accuracy of the simulation can be stabilized at a high level when there is enough training data in the long-period simulation.

Figures A2 and A3 elucidate the spatial distribution characteristics of the state accuracy of the Bayesian network simulation and that of the GREB model simulation, which provides a basis for the subsequent selection of regions for recalculating GREB simulation data based on the Bayesian network simulation result. The state accuracy of the Bayesian network simulation result is relatively uniform in spatial distribution and has no obvious spatial characteristics. However, the state accuracy of the GREB model has obvious characteristics of latitude differentiation. Based on the above, the state accuracy on the space is averaged along the latitudinal direction as shown in Fig. A4. The variances of the state accuracy of the Bayesian network simulation result in six cases are 0.016 (BN-5-T_surf), 0.014 (BN-7-T_surf), 0.014 (BN-9-T_surf), 0.017 (BN-5-T_atmos), 0.008 (BN-7-T_atmos), and 0.004 (BN-9-T_atmos), and the variances of state accuracy of the GREB model simulation result in six cases are 0.089 (GREB-5-T_surf), 0.070 (GREB-7-T_surf), 0.060 (GREB-7-T_surf), 0.077 (GREB-5-T_atmos), 0.054 (GREB-7-T_atmos), and 0.036 (GREB-9-T_atmos). The variance indicates that the fluctuation range of the state accuracy of the Bayesian networks is much smaller than that of the GREB model along the latitude direction. This means that Bayesian networks have a wide range of applications in global climate state simulation.

Although the Bayesian networks have a higher state accuracy in both simulations, we also found that the state simulation accuracy of the GREB model in the range of 30^∘ S to 30^∘ N tends to be higher than that of the Bayesian networks when the classification numbers are 5, 7, and 9 in the surface temperature simulation. Therefore, we think that the GREB model can accurately represent the surface temperature simulation process in this range, and there is no abrupt change in region that cannot be expressed, so in the subsequent optimization, only the range of 30 to 90^∘ N and 30 to 90^∘ S is selected as the optimization region for surface temperature simulation.

Based on the above comparative analysis of the state accuracy of Bayesian network simulations and the state accuracy of GREB model simulations, the range of 30 to 90^∘ N and 30 to 90^∘ S was selected as the empirical parameter for the range of subsequent data recalculating in surface temperature simulations, and the global range was selected as the range of data recalculating in temperature of the atmosphere simulations.

Figure A1 Comparison results of the mean of state accuracy of different methods. The labels 5, 7, and 9 are the different classifications; BN-T_surf represents the mean of state accuracy of surface temperature by the Bayesian network simulation; GREB-T_surf represents the mean of state accuracy of surface temperature by the GREB model; BN-T_atmos represents the mean of state accuracy of temperature of the atmosphere by the Bayesian network simulation; GREB-T_atmos represents the mean of state accuracy of temperature of the atmosphere by the GREB model simulation.

Figure A2 Comparison results of the state accuracy of simulation results of surface temperature. (a) State accuracy of Bayesian networks under five classifications; (b) state accuracy of Bayesian networks under seven classifications; (c) state accuracy of Bayesian networks under nine classifications; (d) state accuracy of the GREB model under five classifications; (e) state accuracy of the GREB model under seven classifications; (f) state accuracy of the GREB model under nine classifications.

Figure A3 Comparison results of the state accuracy of simulation results of temperature of the atmosphere. (a) State accuracy of Bayesian networks under five classifications; (b) state accuracy of Bayesian networks under seven classifications; (c) state accuracy of Bayesian networks under nine classifications; (d) state accuracy of the GREB model under five classifications; (e) state accuracy of the GREB model under seven classifications; (f) state accuracy of the GREB model under nine classifications.

Figure A4 Mean of the spatial distribution of the state accuracy in the latitudinal direction; BN-5 represents the state accuracy of Bayesian networks under the five state classification simulation, BN-7 represents the state accuracy of Bayesian networks under the seven state classification simulation, BN-9 represents the state accuracy of Bayesian networks under the nine state classification simulation, GREB-5 represents the state accuracy of the GREB model under the five state classification simulation, GREB-7 represents the state accuracy of the GREB model under the seven state classification simulation, and GREB-9 represents the state accuracy of the GREB model under the nine state classification simulation. (a) Comparison results of surface temperature; (b) comparison results of temperature of the atmosphere.

Appendix B: Tables

Table B1 Five classifications to indicate climate state.

Table B2 Seven classifications to indicate climate state.

Table B3 Nine classifications to indicate climate state.