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

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.


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

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  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 2 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 simulationtime. 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 32 × 32 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 setups 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 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 RH 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 RH 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, (T SR − OLR in MITgcm atmospheric module, in units of [Wm −2 ]) 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 −2 ]) is estimated in two ways: (i) from the net heat flux into the ocean (T FLUX 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 H per unit area, defined as: where ρ = 1023 kg m −3 is seawater density, c p = 4000 J (K kg) −1 is the specific heat capacity (same values used in Hobbs et al (2016)), T is seawater potential temperature and A 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) . Land mask, vegetation cover, soil albedo and runoff are given to initialise the coupled atmosphereocean model. Simulations are run until deep-ocean equilibrium and the last 100 yr are used for diagnostics.

Results
We describe here the results obtained from the analysis of the suite of modelling set-ups listed in Table 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 RH 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.4 • C 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 (T f in − T init )/T init [ • C] = 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 47 • N, 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 12.27 · 10 6 to 19.85 · 10 6 km 2 in Howell et al (2016)), as can be seen from Fig. 5. However, the Antarctic sea ice cover is too extensive as compared to 12 · 10 6 km 2 in the observations and to the annual average of 20.3 · 10 6 km 2 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.

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 RH which is set to RH = 0.85, 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 −2 , giving rise to SAT that reaches 15.2 • C 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 (T f in − T init )/T init [ • C] = 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 15 · 10 6 km 2 . It turns out that there is a strong dependence on the RH parameter that we have summarised in Fig. 6 by plotting the temperature drift in the ocean as a function of the rel-  ative error in this parameter, (RH − RH c )/RH c , where RH c 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 RH where the drift is practically equal to zero (this value, RH = 0.844, 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 dT /dt 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 2012Gupta et al , 2013. Our analysis shows that choosing a value of RH 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.

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 Run4 Run5 Run1 Fig. 6: Sensitivity to the RH parameter. Drift in the mean ocean temperature vs. the relative error in the RH parameter, (RH − RH c )/RH c , where RH c corresponds to no-drift, for different modelling set-ups: 8-level SPEEDY cloudalbedo parameterisation as for Run1 and Run2 (dots); 5level SPEEDY cloud-albedo parameterisation as for Run3 (squares); friction heating included as in Run4 (diamonds); Callovian simulations as Run5 (stars). 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 RH 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 dT /dt is strongly sensitive to the RH parameter: a huge change in dT /dt (of the same order as that occurring between Run1 and Run2) is obtained for very small variations of RH (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 RH 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 wellstabilised control runs. This criterium, applied to different tuning parameters and/or parameterisations, can be generalised to other coupled climate models.

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 −2 in other coupled atmosphereocean general circulation models, such as HadCM3 and its coarse-resolution version FAMOUS (Pascale et al 2011). Moreover, biases ranging from −0.3 to 4.8 Wm −2 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.8 • C (Fig. 1), and from the global ocean temperature, that decreases of only (T f in − T init )/T init [ • C] = −2% 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 • N (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).

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 re-spect 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 RH 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 presentday 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 2 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 18 • C and 3.8 • C, 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 12 · 10 6 km 2 in the north polar region, while it shows larger variability in the south polar region near a mean value of 9 · 10 6 km 2 (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 RH, 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.

Discussion and conclusions
High-accuracy measurements from satellite programs as CE-RES (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 −2 . The resulting imbalance at TOA is in agreement with that determined from changes in ocean heat content that amounts to 0.8 ± 0.2 Wm −2 (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 −2 .
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 −0.55 Wm −2 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 firstorder 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 T FLUX) 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 • N (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 temper-ature drift dT /dt depends on the tuning parameter RH 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 RH 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 RH 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 atmosphereocean 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.