In this post, we will be discussing the building of a simple climate model. We will then run simple experiments to help us understand the model and to explain how our climate, in general, reacts to certain forcings. In the first section, we will introduce the governing equations and both define and explain each term. Attached in the supplemental information, we will have the code that we built in order to solve the climate evolution. In the second section, we will simulate certain scenarios.
The model
The model that we are using is a two box model. In this model, we have a box for the atmosphere and a box for the surface. Intrinsic in our assumptions, we also have a box for space. We will not model this box, instead we will assume that all the energy comes from the solar input, and exits into space, to never come back. The following model is an unsteady two box model. The variable that we are most interested in is the temperature of the box. Since the model is unsteady, we have two differential equations for the temperature of surface and the atmosphere with respect to time. The equations in general are
\[\begin{array}{*{20}{c}}{\left( 1 \right)}&{\frac{{d{T_{{\rm{surf}}}}}}{{dt}} = {\beta _{{\rm{surf}}}}\left( {\frac{{d{E_{{\rm{surf, in}}}}}}{{dt}} – \frac{{d{E_{{\rm{surf, out}}}}}}{{dt}}} \right)}\end{array}\]
\[\begin{array}{*{20}{c}}{\left( 2 \right)}&{\frac{{d{T_{atms}}}}{{dt}} = {\beta _{{\rm{atms}}}}\left( {\frac{{d{E_{atms{\rm{, in}}}}}}{{dt}} – \frac{{d{E_{atms,{\rm{ out}}}}}}{{dt}}} \right)}\end{array}\]
where T is temperature, β is a constant, E is an energy flux, and t is time. The constant, β, is as follows:
\[{\beta _{{\rm{surf}}}} = \frac{{3600\frac{s}{{hr}} \cdot 24\frac{{hr}}{{day}} \cdot 365.25\frac{{day}}{{yr}}}}{{{\rho _{water}} \cdot {C_{v,water}} \cdot MLD}}\]
\[{\beta _{{\rm{atms}}}} = \alpha \frac{{3600\frac{s}{{hr}} \cdot 24\frac{{hr}}{{day}} \cdot 365.25\frac{{day}}{{yr}}}}{{10,000\frac{{kg}}{{{m^2}}} \cdot {C_{v,air}}}}\]
where ρwater is the density of water, Cv, water is the heat capacity of water, Cv, air is the heat capacity of air, and MLD is the mixed layer depth. In equations (1) and (2), we need to understand the energy fluxes to determine how the temperature evolves on the surface and in the atmosphere. We will begin with the surface.
The surface
The surface energy inflows are modeled with equations (3) and (4). Equation (3) represents the solar influx to the surface.
\[\begin{array}{*{20}{c}}{\left( 3 \right)}&{\frac{{d{E_{{\rm{surf, in 1}}}}}}{{dt}} = \frac{{d{E_{sol}}}}{{dt}}\left( {1 – {R_{{\rm{surf}}}} – {R_{{\rm{atms}}}} – Ab{s_{{\rm{atms}}}}} \right)}\end{array}\]
where dEsol/dt is the solar input, Rsurf is the ratio of energy that is reflected from the surface, Ratms is the ratio of energy that is reflected by the atmosphere, and Absatms is the amount of short wave radiation absorbed by the atmosphere. All parameters are constants in our model except for Rsurf, which is a function of the surface temperature, given by the following equation:
\[{R_{{\rm{surf}}}} = 0.08465 + 0.38{e^{ – 0.006{{\left( {{T_{Esurf}} – 260} \right)}^2}}}\]
Equation (4) represents the amount of longwave radiation that is radiated from the atmosphere.
\[\begin{array}{*{20}{c}}{\left( 4 \right)}&{\frac{{d{E_{{\rm{surf, in 2}}}}}}{{dt}} = 1.25 \cdot em{s_{{\rm{atms}}}} \cdot \sigma \cdot T_{{\rm{atms}}}^4}\end{array}\]
where σ is the Stefan-Boltzmann constant and emsatms is the emission of radiation by the atmosphere. emsatms is given by
\[em{s_{{\rm{atms}}}} = 0.76 + 0.03\sqrt {C{O_2}/320} + 0.1{H_2}O\]
where CO2 is the concentration of carbon dioxide and H2O is the concentration of water. CO2 is a constant in our first models, but H2O is given by the following function:
\[{H_2}O = {K_{{H_2}O}} \cdot \left[ {0.6 + 0.5{e^{5420\left( {{T_{{\rm{surf}}}} – 288.03} \right)/\left( {288.03 \cdot {T_{{\rm{surf}}}}} \right)}}} \right]\]
where KH2O is a constant.
Now, the outflows are modeled with equation (5).
\[\begin{array}{*{20}{c}}{\left( 5 \right)}&{\frac{{d{E_{{\rm{surf, out}}}}}}{{dt}} = \sigma \cdot T_{{\rm{surf}}}^4 + LS{H_{{\rm{flux}}}}}\end{array}\]
where LHSflux is the latent and sensible heat flux moving from the surface to the atmosphere. LHSflux is
\[LS{H_{{\rm{flux}}}} = 104 \cdot \sqrt {\left( {{T_{{\rm{surf}}}} – {T_{{\rm{atms}}}}} \right)/20.55} \]
Here the energy moves from the surface to the atmosphere. The value of LHSflux can only be positive; a negative number to the non-integer exponent is a imaginary number.
The atmosphere
The atmosphere inflow rates are modeled with equations (6) and (7). Equation (6) explains the solar inflow into the atmosphere.
\[\begin{array}{*{20}{c}}{\left( 6 \right)}&{\frac{{d{E_{{\rm{atms, in 1}}}}}}{{dt}} = Ab{s_{{\rm{atms}}}}\frac{{d{E_{sol}}}}{{dt}}}\end{array}\]
Equation (7) models the amount of energy that is radiated from the surface as well as the flux of sensible and latent head from the surface.
\[\begin{array}{*{20}{c}}{\left( 7 \right)}&{\frac{{d{E_{{\rm{atms, in 2}}}}}}{{dt}} = \sigma \cdot em{s_{{\rm{atms}}}} \cdot T_{{\rm{surf}}}^4 + LS{H_{{\rm{flux}}}}}\end{array}\]
Equation (8) represents the emissions from the atmosphere to the surface and space.
\[\begin{array}{*{20}{c}}{\left( 8 \right)}&{\frac{{d{E_{{\rm{atms, out}}}}}}{{dt}} = 2.0 \cdot \sigma \cdot em{s_{{\rm{atms}}}} \cdot T_{{\rm{atms}}}^4}\end{array}\]
The parameters
There are physical parameters for the model as well as computational parameters for the computational portion of the model. We used the following values for our parameters for our base case.
| Parameter | Value | Unit |
|---|---|---|
| ρw | 1000 | kg/m3 |
| Cv,water | 4184 | J/kg K |
| Cv,air | 700 | J/kg K |
| α | 1.48 | - |
| dEsol/dt | 342 | W/m2 |
| Ratms | 0.225 | - |
| Absatms | 0.196 | - |
| CO2 | 320 | ppm |
| σ | 5.67E-8 | W/m2K4 |
| KH2O | 1 | - |
| dt | 0.001 | yrs |
Steady State
Before running the model, initial conditions are necessary for both the temperature of the surface and the atmosphere. Using the initial conditions of 288.99 K for the surface and 267.44 for the atmosphere, we run the model for 30 years. Here the model can reach a steady state equilibrium. The values are close to the initial condition and are 288.07 K and 267.49 K for the surface and atmosphere, respectively. The following figure, (1), shows the evolution.
Test cases
Sensitivity
When we perturb the base parameters, the final steady state solution, i.e. the final temperatures change. In this subsection, we will change the incoming solar input by 3%. Figure (2) shows the model with the same initial conditions as Figure (1) but with a solar input 97% of the base parameter, and Figure (3) with a solar input 103% of the base parameter. In both figures the initial conditions are the steady state equilibrium values seen in the steady state case. As we expected, in Figure (2), we see that the steady state temperature is less than the base level temperatures, and in Figure (3), we see that the steady state temperatures.
If we run multiple tests in between these values, we can see the sensitivity of the models final solution to the solar input. This can be seen in Figure (4). In this figure we see that the slope of the line decreases as the solar input becomes larger. Meaning, the steady state solutions are more sensitive to the solar input
Response time
The model will react to changes in the parameters in a certain time period, response time. As in many unsteady models, the model will reach an equilibrium asymptotically; meaning, it will approach the final state but never reach it (assuming no numerical precision errors). From Figure (2) and (3), we calculated the time it took for the system to reach 95% of the total change. When the solar input is reduced to 97% of its original value it takes 22.7 yrs for the atmosphere to adjust and 22.9 yrs for the surface to adjust. When the solar input was increase to 103% it took the surface 12.6 yrs to adjust and the atmosphere 12.3 yrs to adjust. In Figure (5), we show the response time as a function of the solar input, i.e. the size of the perturbation. We see that as the solar input in larger than the base level, the response time does not change much. Near the base level, there is still a response time. This is not correct, but due to our definition of the response time, it has a larger value. That is, the response time is the amount of time it takes for the model to reach 95% of the difference, regardless of small the difference is. The most import part of this plot is that decrease in solar input shows a longer response time of the model. Meaning, a increase in solar energy will cause the earth to heat rather quickly, but the removal of solar energy will cause the earth to response rather slowly.
Another interesting result we found was about the interaction between the atmosphere and the surface and their respective response times. In Figure (6), we increased the initial temperature of the surface by 10 degrees. Here, we see that the temperatures decrease towards the equilibrium conditions we saw in Figure (1). Interestingly, the temperature of the atmosphere almost instantly adjusts to the surface, even though there was no initial perturbation. This result was also seen in Figure (7). In this figure, the atmosphere was given a +10 degree initial perturbation. One can barely make out the spike in this plot. The atmosphere almost immediate responds and decrease to the equilibrium temperature. This figures tell us that it takes time for the surface to equilibrate to perturbation, and the atmosphere equilibrates almost instantly to the surface. The surface controls the response time, and the atmosphere responds to the surface. Therefore, the atmosphere has nearly the same response time as the surface when the surface is changing.
Cloud cover
In this subsection, we are going to experiment with cloud cover. Here, cloud cover will be a function of the surface temperature, given by the following table. This will add an extra feedback mechanism in our model. For values in between that data points, we will use linear interpolation. From our table, we see that our steady state equilibrium temperature will remain at steady state because the cloud cover is the exact same as our base parameters at this temperature. If we add a perturbation to the system, we can change the steady state and see the effects of the cloud cover feedback. In Figure (8), we added a 103% change to the solar input, just like Figure (3). Compared to Figure (3), we see that Figure (8) shows that the final temperatures are cooler. As the temperature increases, the cloud cover also increases. This causes the amount of solar input to the surface to decrease due to increased reflection from clouds. This buffers the temperature by making the system less sensitive to the increase of solar input.
| surface temperature[K] | cloud cover [-] |
|---|---|
| 258 | 0.015 |
| 268 | 0.11 |
| 278 | 0.36 |
| 288 | 0.6 |
| 298 | 0.84 |
| 308 | 0.92 |
| 318 | 0.95 |
Greenhouse effect
The greenhouse effect allows the earth to reach a warmer temperature at equilibrium. Before we test the effects of the changing our parameters, we will first completely remove the greenhouse effect. Here, we will remove the 1st term in equation (7). This removes the input of the radiation from the surface to the atmosphere, in turn lowering the amount of the radiation of the atmosphere back to the surface. This energy is lost to space, and therefore, the overall temperature of the system, i.e. the earth, will be lower. Figure (9) shows the decrease in the temperatures due to the removal of the greenhouse effect. The temperatures of the earth goes to about 235K and 210K for the surface and atmosphere, respectively.
We can also enhance the effects of the greenhouse effect by increasing the concentration CO2 in our model. Originally, in our base model, we had a CO2 concentration of 320 ppm. By increasing the concentration of CO2, we increase the value of emsatms. When this parameter increases, the amount of energy it radiates back to the surface increases, increasing the greenhouse effect. In Figure (10), we doubled the CO2 concentration to see the effect on the final steady state temperatures. Here, we see that the temperature increases ~2 degrees. This value is with ranges that more sophisticated models predict. Differences in the model results and predictions could be due to more complex feedbacks and mechanisms included in their models or a larger and more resolved domain. Regardless, it is remarkable that a simple model with only two boxes can predict an increase in temperature within a range of the more complex models. However, it can only do just that; other models may be able to to predict changes in weather patter, seasonality, etc. that our model cannot predict.
