How to reduce long-term drift in present-day and deep-time simulations?

How to reduce long-term drift in present-day and deep-time simulations?

Maura Brunetti M. Brunetti GAP-Climate, Institute for Environmental Sciences,
University of Geneva,
66 Bd Carl-Vogt, 1205 Geneva, Switzerland
Tel.: +41-22-379 06 25
Fax: +41-22-379 07 44
22email: maura.brunetti@unige.chC. Vérard Climatic Change and Climate Impacts Group,
Institute for Environmental Sciences, University of Geneva,
66 Bd Carl-Vogt, 1205 Geneva, Switzerland
   Christian Vérard M. Brunetti GAP-Climate, Institute for Environmental Sciences,
University of Geneva,
66 Bd Carl-Vogt, 1205 Geneva, Switzerland
Tel.: +41-22-379 06 25
Fax: +41-22-379 07 44
22email: maura.brunetti@unige.chC. Vérard Climatic Change and Climate Impacts Group,
Institute for Environmental Sciences, University of Geneva,
66 Bd Carl-Vogt, 1205 Geneva, Switzerland
Received: date / Accepted: date

Climate models are often affected by long-term drift that is revealed by the evolution of global variables such as the ocean temperature or the surface air temperature. This spurious trend reduces the fidelity to initial conditions and has a great influence on the equilibrium climate after long simulation times. Useful insight on the nature of the climate drift can be obtained using two global metrics, i.e. the energy imbalance at the top of the atmosphere and at the ocean surface. The former is an indicator of the limitations within a given climate model, at the level of both numerical implementation and physical parameterisations, while the latter is an indicator of the goodness of the tuning procedure. Using the MIT general circulation model, we construct different configurations with various degree of complexity (i.e. different parameterisations for the bulk cloud albedo, inclusion or not of friction heating, different bathymetry configurations) to which we apply the same tuning procedure in order to obtain control runs for fixed external forcing where the climate drift is minimised. We find that the interplay between tuning procedure and different configurations of the same climate model provides crucial information on the stability of the control runs and on the goodness of a given parameterisation. This approach is particularly relevant for constructing good-quality control runs of the geological past where huge uncertainties are found in both initial and boundary conditions. We will focus on robust results that can be generally applied to other climate models.

Tuning Energy budget GCM Paleoclimate
journal: Climate Dynamics

1 Introduction

Paleoclimate simulations of the geological past are particularly challenging since initial conditions are not well constrained by sedimentary data, flawed by uncertainties in dating and spatial scarcity, while boundary conditions are often affected by large uncertainty in paleogeographic reconstructions (National Research Council et al., 2011, p. 86). As a consequence, deep-time simulations are not only limited by the usual biases of climate models, but also by additional biases which come from the construction of imperfect initial and boundary conditions. Moreover, the technique of restoring surface temperature and salinity to observed initial conditions, which is used to improve the stability of coupled simulations in preliminary integrations (Sanchez-Gomez et al, 2016), is not possible when accurate initial conditions are not available. In such a situation the tuning procedure assumes a crucial role (Mauritsen et al, 2012; Golaz et al, 2013; Hourdin et al, 2016). In present-day simulations, the tuning procedure guarantees the construction of control runs, which are runs where the climate forcing (e.g., from solar brightness, atmospheric concentration of greenhouse gases, …) is held constant, with minimal spurious drift (Covey et al, 2006; Sen Gupta et al, 2012, 2013). In paleoclimate simulations, the tuning procedure becomes even more important because it helps in constraining the possible values of parameters, thus reducing the number of climatic attractors that can be explored within a given climate model (Freire et al, 2008).

In general, tuning is not well-documented in climate simulations but the scientific community is now more and more aware that tuning should be made a more explicit process and should be taken into account for evaluating and interpreting model results (Hourdin et al, 2016). Tuning is used to improve the performance of a model in reproducing a given climate. However, if tuning is disconnected from the development of improved physical parameterisations at the process level, the risk is to have heavily tuned models that mask the presence of systematic errors (Jakob, 2014). We will present here a tuning procedure that highlights the goodness of a given parameterisation. The idea is to apply the same tuning procedure to a hierarchy of configurations of a climate model, characterised by a given complexity in the representation of its physical processes. We will take advantage of the modular design of the MIT general circulation model (MITgcm) (Marshall et al, 1997a, b; Adcroft et al, 2004), where, as in many other climate models, physical processes and parameterisations in each component of the climate system are designed as a module that can be activated at runtime, allowing one to change the complexity of the climate simulations by including/excluding different parts of the code. The result is that the link between tuning and process description narrows the size of the parameter space, with clear advantages in particular in the case of paleoclimate simulations.

The tuning procedure can change from one climate code to the other (Gregoire et al, 2011; Irvine et al, 2013; Hourdin et al, 2016). The optimal solution found by tuning only one parameter at a time can differ from the one found by perturbing multiple parameters systematically, using objective methods such as, for example, cost function optimisation (see Hourdin et al (2016) and references therein). Such objective methods are still not commonly implemented in climate groups nowadays, while in general tuning procedure is performed in two stages. In a first stage, the model components (atmosphere, ocean, land, sea ice) are finalised independently in forced mode. In a second stage, the components are coupled together and only few parameters are allowed to change (in order to avoid compensating errors). For example, in CESM (Community Earth System Model) these tuning parameters are the sea ice albedo, and the humidity threshold that controls the formation of low clouds (Gent et al, 2011; Tang et al, 2016). We apply here the same tuning procedure to the case of MITgcm in coupled atmosphere-ocean-sea ice configurations. We consistently use the same procedure at each level of complexity of the model to obtain quasi-equilibrium simulations over long time-scales. We start from considering present-day simulations with different physical parameterisations. This will allow us to better understand the limitations of the code and to set the robustness of the results. Then we move to deep-time paleoclimate simulations in order to set the right procedure for obtaining well-balanced control runs.

In both present-day and deep-time simulations, the tuning parameters are optimised by checking the energy imbalance of the system under consideration at the top of the atmosphere (hereafter named TOA imbalance) and at the ocean surface, since the ocean becomes the dominant energy store in the Earth system on a timescale of about one year (Palmer and McNeall, 2014). We show that useful information on the limitations of the activated modules and on the quality of the control runs can be obtained from these energy budgets.

The paper is organised as follows: in section 2 we describe the coupled simulations and the suite of modelling set-ups. In section 3 we analyse the control runs for each configuration and we describe the method that links tuning and parameterisation. In section 4 we discuss the relevance of the method for deep-time simulations and we draw conclusions.

2 Model description and experiments

The coupled model that we use for the present study is the MIT general circulation model (Marshall et al, 1997a, b; Adcroft et al, 2004). We have chosen this code since it is modular programmed, open-access and it can use cubed-sphere grids that turn out to be particularly useful when polar regions are covered by oceans, as was repeatedly the case in the geological past. In this code, the same dynamical kernel is employed for representing atmosphere and ocean dynamics (Marshall et al, 2004) on the same cubed-sphere grid. The Gent and McWilliams scheme (Gent and McWilliams, 1990) is used to parameterise mesoscale eddies, while the KPP scheme (Large et al, 1994) accounts for vertical mixing processes in the ocean surface boundary layer and the interior. The 5-level SPEEDY physics package for the atmosphere, that is described in Molteni (2003), comprises a four-band radiation scheme, boundary layer and moist convection schemes, resolved baroclinic eddies and diagnostic clouds. Orbital forcing is prescribed at present-day values and the atmospheric CO content is fixed to a constant value of 326 parts per million. The Winton thermodynamic model (Winton, 2000) is used for the sea ice component. A 2-layer land model is also included (Hansen et al, 1983).

The first step in using a climate model is to construct the so-called control run, i.e. the quasi-equilibrium configuration obtained by setting the climate forcing to a constant level. The quasi-equilibrium configuration corresponds to a state where global metrics, such as the surface air temperature (SAT) or the global ocean temperature, reach stationary values so that the climate model does not experience drifts that prevent to study the system over long simulation-time. This quasi-equilibrium configuration is obtained after a spin-up phase. It has been shown that long spin-up has the advantage of reducing the rate of climate drift (Sen Gupta et al, 2012), but on the other hand, even in the presence of a small drift, a long integration means that the climate state is more likely to diverge from the initial conditions (Sen Gupta et al, 2013). From here the importance of constructing control runs with minimal drift from the beginning to obtain at the same time great fidelity to initial conditions and the possibility of performing long-term runs.

Here we consider a low resolution cubed-sphere grid, where each face of the cube has cells, giving a horizontal resolution of ca. 2.8. The ocean has 15 vertical levels with different thickness, from 50 m near the surface to 690 m in the abyss. In this way, we can perform simulations over thousand years, which is the dynamical time-scale of the overturning in the entire ocean, in a reasonable amount of CPU time, allowing us to test for different parameters and configurations. One of the advantages of the MITgcm code is that the modules for a given parameterisation can be activated at runtime so that one can easily construct a hierarchy of models with different degrees of complexity by including or not a given physical process. The simulations considered in the present study according to our suite of modelling set-ups are listed in Table 1.

Run1 is the reference run where, beside all the above listed parameterisations, we also included a modification of the bulk cloud albedo that has been implemented in 8-level SPEEDY (Molteni and Kucharski, 2006; Kucharski et al, 2006) in order to correct a too strong high-latitude solar radiation flux. This reference run has been chosen among a series of numerical experiments where we changed the values of vertical diffusivity within the ocean and of snow/ice albedo (as listed in Table 2) so that the temperature drift in the global ocean temperature was minimal, following the tuning procedure used in other modelling groups (Gent et al, 2011; Tang et al, 2016). These albedo values are within the observed range reported, for example, in Nguyen et al (2011). Afterwards, only the relative humidity threshold for low clouds (a parameter referred as RHCL2 in MITgcm atmospheric module) has been changed for tuning the suite of modelling set-ups and obtaining the corresponding control run.

In Run2 we consider a slight change in the input parameter with respect to the value found for Run1 to show its impact on the energy budget. In Run3 we use a different parameterisation for the bulk cloud albedo (the one typically used in 5-level SPEEDY (Molteni, 2003)) in order to discuss the impact of a different parameterisation on energy budget and tuning procedure, and how to decide which parameterisation is better. In Run4 the kinetic energy dissipated within the ocean and the atmosphere is re-inserted into the system as thermodynamic energy allowing for an improved TOA budget. In Run5 we apply the same procedure used in present-day runs to a different deep-time configuration, namely the Callovian (165 Ma) bathymetry (see Brunetti et al (2015) for an example of Callovian configuration).

Our analysis is based on the calculation of annual means of globally integrated energy budget variables, as suggested in Hobbs et al (2016) for a useful first-order diagnostic of the model behaviour. The TOA energy imbalance is calculated by summing the net shortwave radiation flux and the outgoing longwave radiation flux at each grid point, ( in MITgcm atmospheric module, in units of [Wm]) and then performing area-weighted averages. Positive values correspond to a net radiative warming of the planet. The budget at the ocean surface (in units of [Wm]) is estimated in two ways: (i) from the net heat flux into the ocean ( in MITgcm ocean component, positive if heat enters into the ocean surface and the ocean warms) and (ii) from taking the time derivative of the ocean heat content per unit area, defined as:


where  kg m is seawater density,  J (K kg) is the specific heat capacity (same values used in Hobbs et al (2016)), is seawater potential temperature and is the ocean surface.

The model is spun up from rest, without snow and ice covered oceans. Initial three-dimensional distributions of potential temperature and salinity are derived from the ocean climatological database (annual means) (Levitus et al, 1994; Levitus and Boyer, 1994). Land mask, vegetation cover, soil albedo and runoff are given to initialise the coupled atmosphere-ocean model. Simulations are run until deep-ocean equilibrium and the last 100 yr are used for diagnostics.

Name Bathymetry Description Cloud albedo as in
-level SPEEDY
Run1 Present-day Reference run 0.8440
Run2 Sensitivity to RH 0.8500
Run3 Cloud albedo 0.7677
Run4 Friction heating 0.7900
Run5 Callovian 0.9035
Table 1: List of simulations.

3 Results

We describe here the results obtained from the analysis of the suite of modelling set-ups listed in Table 1.

3.1 Reference run: Run1

The reference run, Run1, is the result of the tuning procedure described in the previous section where the final tuning parameter is the relative humidity threshold for low clouds. Thus, given the approximations of the set-up described in Section 2, this control run can be considered as a good description of the pre-industrial climate.

The SAT reaches an average value of 13.4C during the last 500 yr of simulation, as can be seen in Fig. 1. The global ocean temperature in Run1 shows a very small drift, with an ocean temperature increase of 2% over 1000 yr including the spin-up phase, as shown in Fig. 2. The vertical section of the ocean temperature evolution in Fig. 3a shows that although during the spin-up phase warming occurs near the surface and cooling at depth, implying vertical redistribution of heat from the deep ocean to the surface, afterwards the ocean has constant conditions in all the vertical section, including the deep ocean in the last 300 yr. The overturning circulation is well established in the Atlantic ocean with a maximum intensity of 17.9 Sv at latitude 47N, that are typical values in low-resolution runs (see Fig. 4a).

The Arctic sea ice extent is comparable to annual mean values obtained by other climate models in pre-industrial simulations (ranging from to  km in Howell et al (2016)), as can be seen from Fig. 5. However, the Antarctic sea ice cover is too extensive as compared to  km in the observations and to the annual average of  km obtained in pre-industrial control simulations by CCSM4 (Landrum et al, 2012), probably because the sea ice module in our code does not include dynamical but only thermodynamical effects. However, the reached values are rather constant during the last 500 yr showing that the simulation is stable during this period of time.

Ocean albedo 0.07
Max sea ice albedo 0.64
Min sea ice albedo 0.20
Cold snow albedo 0.85
Warm snow albedo 0.70
Old snow albedo 0.53
Vertical diffusivity s
Table 2: Model parameters used to tune the reference run (Run1) and applied to all the other simulations.
Figure 1: Time evolution of annual-averaged surface air temperature (SAT).
Figure 2: Time evolution of annual-averaged global ocean temperature.
Figure 3: Deviation of the annual-averaged global ocean temperature from the first year as a function of depth.
Figure 4: Atlantic meridional overturning circulation (AMOC) in Sv  m/s as a function of depth and latitude, averaged over the last 100 yr of simulation time. (a) Run1. AMOC anomalies with respect to Run1 for the simulations (b) Run2, (c) Run3 and (d) Run4.
Figure 5: Time evolution of annual-averaged sea ice extent.

3.2 Run2

Run2 is analysed here to illustrate the impact of tuning on the simulation results. The same input parameters used in Run1 are employed for this simulation except for the value of which is set to , a value only 0.7% larger than the one used in Run1 (see Table 1). The effect of this small increase is that the energy budget at the ocean surface is unbalanced, 0.25 Wm, giving rise to SAT that reaches 15.2C at the end of the simulation, as shown in Fig. 1, and to a huge drift in the global ocean temperature, with an increase of 23% over 1000 yr (see Fig. 2). In the vertical section (Fig. 3b) the ocean temperature increases at all depths (and especially within the first 1000 m depth), while both the intensity of the Atlantic Meridional Overturning Circulation (AMOC) and the sea ice extent decrease with respect to the reference run, as shown in Figs. 4b and 5, respectively. We remark that the Antarctic sea ice, continuously decreasing throughout the simulation, is more sensitive to drift than the Arctic sea ice, which stabilises to an equilibrium value of  km.

It turns out that there is a strong dependence on the parameter that we have summarised in Fig. 6 by plotting the temperature drift in the ocean as a function of the relative error in this parameter, , where corresponds to the no-drift case. Dots in this plot refer to runs with the cloud-albedo parameterisation as in 8-level SPEEDY. We can see that there is a value of where the drift is practically equal to zero (this value, , corresponds to the reference run, Run1). Small variations around this value gives rise to unbalanced runs, like Run2.

Note that drift in the global ocean temperature occurs in many of the state-of-the-art climate model simulations conducted under CMIP5 (Hobbs et al, 2016) where is typically higher than the value reached in Run1. Given that models drift away from the observed state chosen at initialisation if tuning is not sufficiently accurate, short spin-up results in greater fidelity to initial conditions with respect to long spin-up. On the other hand, long spin-up is an important factor in drift reduction (Sen Gupta et al, 2012, 2013). Our analysis shows that choosing a value of that minimises the drift, such in Run1, increases model fidelity to initial conditions and at the same time guarantees stability over simulation times of the order of thousand years.

Simulation TOA imbalance Surface imbalance
[Wm] [Wm] [K/Century]
Run1 2.65 0.00 0.001
Run2 2.81 0.25 0.052
Run3 2.74 0.23 0.047
Run4 -0.55 0.04 0.009
Run5 -0.36 0.09 0.016
Table 3: Energy imbalances and drift in the global ocean temperature.
Figure 6: Sensitivity to the parameter. Drift in the mean ocean temperature vs. the relative error in the parameter, , where corresponds to no-drift, for different modelling set-ups: 8-level SPEEDY cloud-albedo parameterisation as for Run1 and Run2 (dots); 5-level SPEEDY cloud-albedo parameterisation as for Run3 (squares); friction heating included as in Run4 (diamonds); Callovian simulations as Run5 (stars).

3.3 Run3

We have tested in Run3 a different way of representing the bulk cloud albedo, to be compared to the one used in Run1. As in the atmospheric module 5-level SPEEDY (Molteni, 2003), the bulk cloud albedo is held constant in Run3. Since this parameterisation gives rise to a solar radiation that is too strong at high latitudes in coupled models, developers of the SPEEDY code remarked that using a bulk cloud albedo that increases with latitude improved the simulation results (Kucharski et al, 2016). We thus expect that the parameterisation used in Run3 is less accurate than the one used in the reference run (Run1).

We have checked that the parameterisation used in Run1 produces less net solar radiation at high latitudes (see Fig. 7) with respect to Run3, as expected from the different description of bulk cloud albedo employed in 8-level and 5-level SPEEDY. We find that while SAT is rather constant during the first 800 yr in Run3, it tends to increase in the last century (see Fig. 1). The annual-averaged global ocean temperature has a positive trend all along the simulation (see Fig. 2). This increase in temperature can also be observed in Fig. 3c where the ocean temperature is seen to raise, especially in the upper 500 m. This change in the ocean and surface air temperature can be related to a decrease of sea ice extent in both north and south polar regions with respect to the reference run, as can be seen in Fig. 5. The AMOC is only slightly affected, as can be seen in Fig. 4c, while the total heat transport is smaller than in the case of Run1, as can be seen in Figs. 8a.

The monotonic increase in global ocean temperature demonstrates that Run3 is not well stabilised. It is important to note that if the tuning parameter is slightly reduced, the result is a monotonic decrease in the global ocean temperature, showing that it is very difficult to obtain a stable control run for this series of simulations. In order to better understand this point, we have investigated the behaviour of the temperature drift as a function of the relative error in parameter (Fig. 6, squares) and interestingly we found that the dependence on this relative error is different from the case of Run1-Run2 series (Fig. 6, dots). In the present case, the temperature drift is strongly sensitive to the parameter: a huge change in (of the same order as that occurring between Run1 and Run2) is obtained for very small variations of (of the order of 0.01%, to be compared with 0.7% in the case of Run1 and Run2). This means that the parameterisation considered in Run3 and all the runs in the same series represented by squares in Fig. 6 is much more sensitive to the tuning parameter than the parameterisation used in Run1 series (dots in Fig. 6). This sensitivity can be used as a criterium to establish for the goodness of a given parameterisation since the more the sensitivity to the tuning parameter, the smaller the range where the tuning parameter can vary to obtain minimal drift and well-stabilised control runs. This criterium, applied to different tuning parameters and/or parameterisations, can be generalised to other coupled climate models.

Figure 7: Zonal mean of TOA net shortwave radiation in simulations Run1 and Run3.

3.4 Run4

The origin of the imbalance at TOA has received a great deal of attention from the scientific community (Pascale et al, 2011; Lucarini and Ragone, 2011; Hobbs et al, 2016). The imbalance is ascribed to physical processes that have been neglected or approximated in climate models. In our set-up, the main source/sink of heat are (i) the neglect of the heating due to kinetic energy dissipation by internal stress and viscous processes (Lucarini and Pascale, 2014) ; (ii) the approximation of considering fixed soil moisture; (iii) the neglect of the sea ice dynamics; (iv) the approximation used in the thermodynamical module to limit the ice thickness to a maximal value of 10 m. Other sources of imbalance have numerical origin and can be related to the time-stepping method, to the advection schemes or to hyperdiffusion operators introduced in the dynamical core in order to smooth and avoid divergences (Lucarini and Pascale, 2014).

The heating due to kinetic energy dissipation is in general dominant with respect to the other effects and can reach values of the order of 3 Wm in other coupled atmosphere-ocean general circulation models, such as HadCM3 and its coarse-resolution version FAMOUS (Pascale et al, 2011). Moreover, biases ranging from to 4.8 Wm are found in the net global balance of pre-industrial simulations included in the IPCC4AR (Lucarini and Ragone, 2011). In order to better understand this point, we have performed simulations where the friction heating is returned back to the atmospheric component by activating the corresponding switch in MITgcm. We obtain a simulation, Run4, where the TOA imbalance is much reduced with respect to the reference run, as can be seen in Table 3. The tuning procedure gives rise to a very well equilibrated simulation, as can be inferred from the time evolution of SAT, that reaches a rather constant value of 12.8C (Fig. 1), and from the global ocean temperature, that decreases of only % in 1100 yr (Fig. 2). From the vertical section of the ocean temperature shown in Fig. 3d, we observe that there is less warming of the upper ocean than in Run1, the AMOC increases almost everywhere by 4 Sv (Fig.  4d) while the sea ice extent becomes very well stabilised in the last 500 yr of the simulation (Fig 5).

The range where the tuning parameter gives reasonable well equilibrated runs is comparable to Run1 series (compare the slope of the curves with dots and diamonds in Fig. 6), thus we can conclude that including friction heating back into the atmospheric column has a global positive effect on the quality of the simulations, since the TOA imbalance is reduced with respect to Run1, and this is a result that can be generalised to other coupled climate models. The positive effect of accounting for the friction heating clearly appears in the total heat transport, that is much closer to observed values of the order of 6 PW at 30(Fasullo and Trenberth, 2008) in Run4 than in all the other considered simulations, as shown in Fig. 8. The improvement is particularly effective in the atmospheric contribution to the total heat transport especially within the northern hemisphere, while it is negligible in the oceanic contribution (compare Run1 and Run4 in Figs. 8a and b).

Figure 8: (a) Total heat transport in PW ( W) and the oceanic contribution to this total transport. A positive value on the vertical axis corresponds to a northward transport. (b) Zoom of the oceanic contribution.

3.5 Run5

We now apply the previous method to the case of a different paleogeographic configuration. Since we have seen that including friction heating back into the atmospheric component had positive effects on the simulation results, we employ the same procedure here. We also use the cloud albedo parameterisation of Run1 that gives improvements with respect to the one used in Run3. The paleogeographic configuration that we consider here corresponds to the Callovian-Oxfordian transition (165 Ma) (see Fig. 9a or, for more details, Fig. 1 in Brunetti et al (2015)), where Pangea breakup gave rise to the formation of the Central Atlantic and the Proto-Caribbean basin, and provided a connection between the Neo-Tethys and Panthalassa. We have already used this configuration in ocean-only simulations presented in Brunetti et al (2015). Apart from the interest of studying this period using a coupled climate model and comparisons with geological data, that we will pursue in forthcoming publications, the technical challenge is now to obtain a stable simulation by tuning the parameter so that the TOA energy imbalance is minimal and the surface energy imbalance is nearly zero. For the present test we do not put much effort in optimising the initial conditions, for which we took the zonal average of the potential temperature over the present-day Pacific ocean, a constant salinity distribution at 39 psu, a vegetation cover that is homogeneous in latitude from Sellwood and Valdes (2008), and a runoff map based on the topography used in Brunetti et al (2015). The initial conditions are anyway affected by many uncertainties, the range of published estimates of Jurassic atmospheric CO content varying between present-day values to 15 times such values (Sellwood and Valdes, 2008).

The simulation results are shown in Fig. 9. Both SAT and global ocean temperature increase during the spin-up phase and reach a stable value around 18C and 3.8C, respectively (Fig. 9b). The vertical section of the ocean temperature shows that it is well stabilised at all depths (Fig. 9c). The upper ocean is much warmer than at the beginning, showing that the initial conditions should indeed be optimised in order to reduce the gap with respect to the final equilibrium. The sea ice extent reaches a constant value of  km in the north polar region, while it shows larger variability in the south polar region near a mean value of  km (Fig. 9d).

The sensitivity to the tuning parameter, quantified by the slope in the plot in Fig. 6 showing temperature drift against the relative error in , is of the same order as the sensitivity in the Run4 series (compare stars and diamonds in Fig. 6). We have thus proved that the above method is general and can be applied successfully to other bathymetric configurations.

Figure 9: Callovian simulation Run5. (a) paleogeography. (b) Time evolution of annual-averaged surface air temperature and global ocean temperature. (c) Deviation of the annual-averaged global ocean temperature from the first year as a function of depth. (d) Time evolution of annual-averaged sea ice extent.

4 Discussion and conclusions

High-accuracy measurements from satellite programs as CERES (Clouds and the Earth’s Radiant Energy System (Wielicki et al, 1996)) and SORCE (Solar Radiation and Climate Experiment (Anderson and Cahalan, 2005)) have constrained the radiation fluxes at TOA with uncertainty range of less than 1 Wm. The resulting imbalance at TOA is in agreement with that determined from changes in ocean heat content that amounts to  Wm (Hansen et al, 2011; Wild et al, 2013). Thus, in the present-day climate, observations point to equal values of TOA and ocean surface budget of the order of 1 Wm.

Ideally the same imbalance should be obtained in climate models. However, while ensemble averages are in general in agreement with these values (Wild et al, 2013; Trenberth et al, 2014), there is a huge spread of results in each model of the ensemble, even in pre-industrial control simulations (see, for example, Fig. 2 in Lucarini and Ragone (2011) for CMIP3 and Fig. 1 in Hobbs et al (2016) for CMIP5). Here, the imbalance should ideally be zero (Lucarini and Ragone, 2011) since these simulations represent an estimate of the unforced quasi-equilibrium climate that evolves under the only effect of internal nonlinear dynamics. Moreover, the values of imbalance at TOA and at ocean surface are in general different in climate models (Hobbs et al, 2016).

In the control runs presented here, the ocean surface budget is tuned to zero to avoid temperature drift. Therefore the imbalance at the ocean surface can be considered as a measure of the goodness of the tuning procedure in control runs. In our simulations, the TOA budget is different from zero, as shown in Table 3. The TOA imbalance can be reduced to lower values by improving the climate code, as in our case where it moved from 2.65 to  Wm when friction heating was taken into account in the simulations (compare Run1 and Run4 in Table 3). Therefore, the TOA budget can be considered as a measure of the goodness of parameterisations and coding, since it is related to the presence of energy sources/sinks within the coupled climate system due to imperfect representations of physical processes (such as, for example, friction heating) and/or to numerical diffusion. Thus, we can confirm the importance of using global metrics, such as TOA and ocean surface imbalance, as first-order diagnostics of the model behaviour (Hobbs et al, 2016). These two quantities should always be explicitly stated when numerical results are presented.

A model can perfectly conserve energy between its components at the ocean surface, as the MITgcm (Campin et al, 2008), and still have energy sources/sinks that affect the TOA budget (due to approximate processes within the model, like numerical or kinetic-energy dissipation), that are not accounted for and are sometimes called ‘ghost energy’ (Lucarini and Ragone, 2011). Within certain models, a correction is applied to eliminate the kinetic energy dissipation by the internal stresses at each time step (Pascale et al, 2011), while in others, analysis of the results are presented as anomalies with respect to the time-mean nonconservation term (Hobbs et al, 2016). We have chosen to explicitly show the values of TOA imbalance since we think that they are very useful to gain insight into the limitations/biases of a given climate model.

Since we have verified that the net ocean surface heat flux (that is an output of the MITgcm called ) exactly accounts for the temperature drift in the ocean volume, we can state that the energy leaks in MITgcm are mainly outside the ocean, as occurs in the majority of climate models considered in Hobbs et al (2016). This is also confirmed by the fact that the oceanic contribution to the total heat transport, shown in Fig. 8b, does not change much for different modelling set-ups, while the total heat transport becomes more closer to the observed value of 6 PW at 30(Fasullo and Trenberth, 2008) when the correct set-up is implemented (compare Run1 to Run4 in Fig. 8a).

We have shown that, taken together, the tuning procedure and the parameterisations can give useful insight into limitations of a climate model. An example for such an approach is given by Fig. 6 where we plotted how the temperature drift depends on the tuning parameter for different configurations of MITgcm. A strong dependence of the temperature drift on the tuning parameter corresponds to a steep growth (as occurs in Run3 series, squares), a small change in giving rise to a huge temperature drift. This means that the implemented physical process strongly depends on the tuning procedure, with the risk that the resulting control run will be not stable and too strongly dependent on internal variability. On the contrary, if small changes in produce the same nearly-zero values for the temperature drift, we can expect to obtain good-quality control runs with the considered set-up of the climate code. Thus, the analysis of the dependence of the temperature drift on the tuning parameter provides insight on the quality of control runs and their stability. This is a general and robust result that can be applied to other climate codes and tuning parameters.

The importance of developing simple diagnostics to assure the quality of control runs and model verification is even more important in simulations of the deep-past where huge uncertainties exist in initial and boundary conditions. From the analysis performed in the previous section, it turns out that the tuning procedure should be applied for each different paleogeographic configuration. We stress here the importance of explicitly stating the values of the TOA and ocean surface energy imbalance, and how the latter depends on the tuning parameter, each time a new simulation of the deep-past is presented. This procedure will indeed allow other researchers to know important aspects on the quality of the numerical results and on the stability of the control run. Unfortunately, nowadays these global metrics are almost never mentioned in paleoclimate studies. We have seen that a small amount of imbalance in the ocean surface budget can induce large effects in SAT, the global ocean temperature or the intensity of the overturning circulation (compare Run1 with Run2 in Figs. 1, 2 and 4). This is particularly important in paleoclimate simulations that need to be run for thousands of years to attain a fully equilibrated coupled atmosphere-ocean state that is not known a priori (while in present-day simulations we know the characteristics that a control run should satisfy at the end of the simulation). Even a small imbalance, but lasting for a long period of time, can strongly affect the final results. This is why it is important in paleoclimate simulations to correctly tune the parameters from the beginning and use procedures, such as the one illustrated in Fig. 6, that allow one to optimise the physical parameterisations and configurations used in the climate code.

In climate simulations of the deep-past there is of course also the need of constructing realistic initial and boundary conditions that can reduce the size of the parameter space. Interdisciplinary effort, with contributions from geologists, physicists, climate modellers, is essential in this respect and intercomparison projects, such as the Paleoclimate Modeling Intercomparison Project (Braconnot et al, 2011; Lunt et al, 2017), are crucial to progress in this field. Sensitivity studies to the initialisation and to the paleogeography are strongly encouraged. Inclusion of the main variables used in the present study (in order to calculate TOA and ocean surface energy imbalance) are recommended in intercomparison datasets.

The computations were performed at University of Geneva on the Baobab and CLIMDAL3 clusters. We thank Jean-Michel Campin, Marjorie Perroud and Martin Beniston for useful discussions, and the MITgcm-support mailing list for valuable advice on the code. This work was partly supported by CTI 15574.1 PFES-ES.


  • Adcroft et al (2004) Adcroft A, Campin JM, Hill C, Marshall J (2004) Implementation of an Atmosphere Ocean General Circulation Model on the Expanded Spherical Cube. Monthly Weather Review 132:2845, DOI 10.1175/MWR2823.1
  • Anderson and Cahalan (2005) Anderson DE, Cahalan RF (2005) The Solar Radiation and Climate Experiment (SORCE) mission for the NASA Earth Observing System (EOS). Solar Physics 230:3–6
  • Braconnot et al (2011) Braconnot P, Harrison SP, Otto-Bliesner B, Abe-Ouchi A, Jungclaus J, Peterschmitt JY (2011) The Paleoclimate Modeling Intercomparison Project contribution to CMIP5. CLIVAR Exchanges No 56 16:15–19
  • Brunetti et al (2015) Brunetti M, Vérard C, Baumgartner PO (2015) Modeling the Middle Jurassic ocean circulation. Journal of Palaeogeography 4:373–386
  • Campin et al (2008) Campin JM, Marshall J, Ferreira D (2008) Sea ice-ocean coupling using a rescaled vertical coordinate z. Ocean Modelling 24:1–14, DOI 10.1016/j.ocemod.2008.05.005
  • Covey et al (2006) Covey C, Gleckler PJ, Phillips TJ, Bader DC (2006) Secular trends and climate drift in coupled ocean-atmosphere general circulation models. J Geophys Res 111:D03,107, DOI 10.1029/2005JD006009
  • Fasullo and Trenberth (2008) Fasullo JT, Trenberth KE (2008) The Annual Cycle of the Energy Budget. Part II: Meridional Structures and Poleward Transports. Journal of Climate 21:2313, DOI 10.1175/2007JCLI1936.1
  • Freire et al (2008) Freire JG, Bonatto C, DaCamara CC, Gallas JAC (2008) Multistability, phase diagrams, and intransitivity in the Lorenz-84 low-order atmospheric circulation model. Chaos 18(3):033121, DOI 10.1063/1.2953589
  • Gent and McWilliams (1990) Gent PR, McWilliams JC (1990) Isopycnal Mixing in Ocean Circulation Models. Journal of Physical Oceanography 20:150–160
  • Gent et al (2011) Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, Lawrence DM, Neale RB, Rasch PJ, Vertenstein M, Worley PH, Yang ZL, Zhang M (2011) The Community Climate System Model Version 4. Journal of Climate 24:4973–4991, DOI 10.1175/2011JCLI4083.1
  • Golaz et al (2013) Golaz JC, Golaz JC, Levy H (2013) Cloud tuning in a coupled climate model: Impact on 20th century warming. Geophysical Research Letters 40:2246–2251, DOI 10.1002/grl.50232
  • Gregoire et al (2011) Gregoire LJ, Valdes PJ, Payne AJ, Kahana R (2011) Optimal tuning of a gcm using modern and glacial constraints. Climate Dynamics 37(3):705–719, DOI 10.1007/s00382-010-0934-8
  • Hansen et al (1983) Hansen J, Russell G, Rind D, Stone P, Lacis A, Lebedeff S, Ruedy R, Travis L (1983) Efficient three-dimensional global models for climate studies: model I and II. Monthly Weather Review 111:609–662
  • Hansen et al (2011) Hansen J, Sato M, Kharecha P, von Schuckmann K (2011) Earth’s energy imbalance and implications. Atmospheric Chemistry & Physics 11:13,421–13,449, DOI 10.5194/acp-11-13421-2011
  • Hobbs et al (2016) Hobbs W, Palmer MD, Monselesan D (2016) An Energy Conservation Analysis of Ocean Drift in the CMIP5 Global Coupled Models*. Journal of Climate 29:1639–1653, DOI 10.1175/JCLI-D-15-0477.1
  • Hourdin et al (2016) Hourdin F, Mauritsen T, Gettelman A, Golaz J, Balaji V, Duan Q, Folini D, Ji D, Klocke D, Qian Y, Rauser F, Rio C, Tomassini L, Watanabe M, Williamson D (2016) The art and science of climate model tuning. Bull Am Meteorological Soc 97:0, DOI 10.1175/BAMS-D-15-00135.1
  • Howell et al (2016) Howell FW, Haywood AM, Otto-Bliesner BL, Bragg F, Chan WL, Chandler MA, Contoux C, Kamae Y, Abe-Ouchi A, Rosenbloom NA, Stepanek C, Zhang Z (2016) Arctic sea ice simulation in the PlioMIP ensemble. Climate of the Past 12:749–767, DOI 10.5194/cp-12-749-2016
  • Irvine et al (2013) Irvine PJ, Gregoire LJ, Lunt DJ, Valdes PJ (2013) An efficient method to generate a perturbed parameter ensemble of a fully coupled aogcm without flux-adjustment. Geoscientific Model Development 6(5):1447–1462, DOI 10.5194/gmd-6-1447-2013, URL
  • Jakob (2014) Jakob C (2014) Going back to basis. Nature Climate Change 4:1042–1045
  • Kucharski et al (2016) Kucharski F, Ikram F, Molteni F, Farneti R, Kang IS, No HH, King MP, Giuliani G, Mogensen K (2016) Atlantic forcing of Pacific decadal variability. Climate Dynamics 46:2337–2351, DOI 10.1007/s00382-015-2705-z
  • Kucharski et al (2006) Kucharski F, Molteni F, Bracco A (2006) Decadal interactions between the western tropical Pacific and the North Atlantic Oscillation. Climate Dynamics 26:79–91, DOI 10.1007/s00382-005-0085-5
  • Landrum et al (2012) Landrum L, Holland MM, Schneider DP, Hunke E (2012) Antarctic Sea Ice Climatology, Variability, and Late Twentieth-Century Change in CCSM4. Journal of Climate 25:4817–4838, DOI 10.1175/JCLI-D-11-00289.1
  • Large et al (1994) Large WG, McWilliams JC, Doney SC (1994) Oceanic vertical mixing: a review and a model with a nonlocal boundary layer parameterization. Reviews of Geophysics 32:363–404, DOI 10.1029/94RG01872
  • Levitus and Boyer (1994) Levitus S, Boyer T (1994) World Ocean Atlas 1994. Volume 4. Temperature. U.S. Department of Commerce, Washington, DC
  • Levitus et al (1994) Levitus S, Burgett R, Boyer T (1994) World Ocean Atlas 1994. Volume 3. Salinity. U.S. Department of Commerce, Washington, DC
  • Lucarini and Pascale (2014) Lucarini V, Pascale S (2014) Entropy production and coarse graining of the climate fields in a general circulation model. Climate Dynamics 43:981–1000, DOI 10.1007/s00382-014-2052-5, 1304.3945
  • Lucarini and Ragone (2011) Lucarini V, Ragone F (2011) Energetics of Climate Models: Net Energy Balance and Meridional Enthalpy Transport. Reviews of Geophysics 49:RG1001, DOI 10.1029/2009RG000323
  • Lunt et al (2017) Lunt DJ, Huber M, Anagnostou E, Baatsen MLJ, Caballero R, DeConto R, Dijkstra HA, Donnadieu Y, Evans D, Feng R, Foster GL, Gasson E, von der Heydt AS, Hollis CJ, Inglis GN, Jones SM, Kiehl J, Kirtland Turner S, Korty RL, Kozdon R, Krishnan S, Ladant JB, Langebroek P, Lear CH, LeGrande AN, Littler K, Markwick P, Otto-Bliesner B, Pearson P, Poulsen CJ, Salzmann U, Shields C, Snell K, Stärz M, Super J, Tabor C, Tierney JE, Tourte GJL, Tripati A, Upchurch GR, Wade BS, Wing SL, Winguth AME, Wright NM, Zachos JC, Zeebe RE (2017) The DeepMIP contribution to PMIP4: experimental design for model simulations of the EECO, PETM, and pre-PETM (version 1.0). Geoscientific Model Development 10:889–901, DOI 10.5194/gmd-10-889-2017
  • Marshall et al (1997a) Marshall J, Adcroft A, Hill C, Perelman L, Heisey C (1997a) A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers. Journal of Geophysical Research 102:5753–5766, DOI 10.1029/96JC02775
  • Marshall et al (1997b) Marshall J, Hill C, Perelman L, Adcroft A (1997b) Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. Journal of Geophysical Research 102:5733–5752, DOI 10.1029/96JC02776
  • Marshall et al (2004) Marshall J, Adcroft A, Campin JM, Hill C (2004) Atmosphere-ocean modeling exploiting fluid isomorphisms. Monthly Weather Review 132:2882–2894
  • Mauritsen et al (2012) Mauritsen T, Stevens B, Roeckner E, Crueger T, Esch M, Giorgetta M, Haak H, Jungclaus J, Klocke D, Matei D, Mikolajewicz U, Notz D, Pincus R, Schmidt H, Tomassini L (2012) Tuning the climate of a global model. Journal of Advances in Modeling Earth Systems 4:M00A01, DOI 10.1029/2012MS000154
  • Molteni (2003) Molteni F (2003) Atmospheric simulations using a GCM with simplified physical parametrizations. I: Model climatology and variability in multidecadal experiments. Climate Dynamics 20:175–191
  • Molteni and Kucharski (2006) Molteni F, Kucharski F (2006) Description of the ictp agcm (speedy), version 40. Tech. rep., Abdus Salam International Centre for Theoretical Physics, Trieste, Italy, URL
  • National Research Council et al. (2011) National Research Council et al (2011) Understanding Earth’s Deep Past: Lessons for Our Climate Future. The National Academies Press
  • Nguyen et al (2011) Nguyen AT, Menemenlis D, Kwok R (2011) Arctic ice-ocean simulation with optimized model parameters: Approach and assessment. Journal of Geophysical Research (Oceans) 116:C04025, DOI 10.1029/2010JC006573
  • Palmer and McNeall (2014) Palmer MD, McNeall DJ (2014) Internal variability of Earth’s energy budget simulated by CMIP5 climate models. Environ Res Lett 9:034,016, DOI 10.1088/1748-9326/9/034016
  • Pascale et al (2011) Pascale S, Gregory JM, Ambaum M, Tailleux R (2011) Climate entropy budget of the HadCM3 atmosphere-ocean general circulation model and of FAMOUS, its low-resolution version. Climate Dynamics 36:1189–1206, DOI 10.1007/s00382-009-0718-1
  • Sanchez-Gomez et al (2016) Sanchez-Gomez E, Cassou C, Ruprich-Robert Y, E F, Terray L (2016) Drift dynamics in a coupled model initialized for decadal forecasts. Climate Dynamics 46:1819–1840
  • Sellwood and Valdes (2008) Sellwood BW, Valdes PJ (2008) Jurassic climates. Proceedings of the Geologists’ Association 119:5–17
  • Sen Gupta et al (2012) Sen Gupta A, Muir LC, Brown JN, J PS, Durack PJ, Monselesan D, Wijffels SE (2012) Climate drift in the CMIP3 models. J Climate 25:4621–4640
  • Sen Gupta et al (2013) Sen Gupta A, Jourdain NC, Brown JN, Monselesan D (2013) Climate drift in the CMIP5 models. J Climate 26:8597–8615
  • Tang et al (2016) Tang Y, Li L, Dong W, Wang B (2016) Reducing the climate shift in a new coupled model. Sci Bull 61:488–494, DOI 10.1007/s11434-016-1033-y
  • Trenberth et al (2014) Trenberth KE, Fasullo JT, Balmaseda MA (2014) Earth’s Energy Imbalance. Journal of Climate 27:3129–3144, DOI 10.1175/JCLI-D-13-00294.1
  • Wielicki et al (1996) Wielicki BA, Barkstrom BR, Harrison EF, Lee RB III, Smith GL, Cooper JE (1996) Clouds and the Earth’s Radiant Energy System (CERES): An Earth Observing System Experiment. Bulletin of the American Meteorological Society 77:853–868
  • Wild et al (2013) Wild M, Folini D, Schär C, Loeb N, Dutton EG, König-Langlo G (2013) The global energy balance from a surface perspective. Climate Dynamics 40:3107–3134, DOI 10.1007/s00382-012-1569-8
  • Winton (2000) Winton M (2000) A Reformulated Three-Layer Sea Ice Model. Journal of Atmospheric and Oceanic Technology 17:525–531
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description