Evolving CO2 Rather Than SST Leads to a Factor of Ten Decrease in GCM Convergence Time

Abstract The high computational cost of Global Climate Models (GCMs) is a problem that limits their use in many areas. Recently an inverse climate modeling (InvCM) method, which fixes the global mean sea surface temperature (SST) and evolves the CO2 mixing ratio to equilibrate climate, has been implemented in a cloud‐resolving model. In this article, we apply InvCM to ExoCAM GCM aquaplanet simulations, allowing the SST pattern to evolve while maintaining a fixed global‐mean SST. We find that InvCM produces the same climate as normal slab‐ocean simulations but converges an order of magnitude faster. We then use InvCM to calculate the equilibrium CO2 for SSTs ranging from 290 to 340 K at 1 K intervals and reproduce the large increase in climate sensitivity at an SST of about 315 K at much higher temperature resolution. The speedup provided by InvCM could be used to equilibrate GCMs at higher spatial resolution or to perform broader parameter space exploration in order to gain new insight into the climate system. Additionally, InvCM could be used to find unstable and hidden climate states, and to find climate states close to bifurcations such as the runaway greenhouse transition.


Groups Runs Descriptions
InvCM-ClimateSensitivity 51 The global mean SST is fixed to 290, 291,  E , 340 K and the 2 CO E mixing ratio evolves to equilibrate the system (see Section2 for details). To obtain the initial SST pattern and atmospheric conditions, we heat the slab ocean and the atmosphere of a modern climate uniformly by the difference between the fixed global-mean temperature and that of the modern climate. The initial 2 CO E mixing ratio is  1 18995 gg E , which is the equilibrium 2 CO E mixing ratio of a preliminary InvCM simulation under 310 K that is not listed in this table. Here, the modern climate means a north-south symmetric and zonally uniform climate state, the SST pattern of which is lowest near the poles (   20.3 C E ), is highest on the equator (  26.3 C E ), and has a global mean value of  11.9 C E . The 2 CO E timescale is 240 days and the mixed layer depth is 50 m. Each simulation runs for 15 yr. The equilibrium 2 CO E as a function of the surface temperature is calculated via a polynomial regression of the geometric mean 2 CO E of the last 10 yr with the random errors smoothed (See Section3.2 and Figure 6). Then we continue to run these simulations for 20 yr under the equilibrium 2 CO E for the corresponding surface temperatures to examine whether equilibrium has been reached.

NormalSlabOcean
3 As in conventional slab-ocean simulations: the 2 CO E mixing ratio is fixed and the atmosphere is coupled to a slab ocean with a uniform depth of 50 m. We choose the 2 CO E mixing ratios to be 18,407, 26,212, and 45,241  1 gg E based on the results from the InvCM-ClimateSensitivity group so that the equilibrium global-mean SST should be approximately 310, 320, and 330 K. The initial condition is a modern climate with a global surface temperature of 289 K. In order to be consistent with other simulations, the thermodynamic properties of the dry air are also taken as those of 2 N E like in InvCM.

InvCM-Ensemble 10
The same equilibrium technique is used as the InvCM-ClimateSensitivity group. The global-mean SST is set to 320 K for all simulations, and we modify the initial SST pattern for the 320 K simulation in the InvCM-ClimateSensitivity group by adding 0 1 . cos cos( , different for each case, so that this group constitutes an ensemble. We run each of the simulations for 10 yr.
The same equilibrium technique is used as the InvCM-ClimateSensitivity group. The mixed layer depth is 50.0 m and the 2 CO E timescale is 30, 60, 120, and 240 days. The global mean SST is set to 310 K for these 4 simulations. For the initial condition, we use the equilibrated 310 K simulation in the InvCM-ClimateSensitivity group.

InvCM-MixedLayerDepth 4
The same equilibrium technique is used as the InvCM-ClimateSensitivity group. The 2 CO E timescale is 240 days and the mixed layer depth is 0.5, 5, 50 and 100 m. The global mean SST is set to 310 K for these 4 simulations. We use the same initial condition as the 310 K simulation in the InvCM-ClimateSensitivity group.

InvCM-InitCondSSTEquilibrium 2
The same equilibrium technique is used as the InvCM-ClimateSensitivity group. The 2 CO E timescale is 240 days and the mixed layer depth is 50 m. The global mean SST is set to 310 K. One simulation (referred to as the control simulation) uses the same initial condition as the 310 K simulation in the InvCM-ClimateSensitivity group. Another simulation (referred to as the experimental simulation) starts from the equilibrium SST pattern obtained from the control simulation and a 2 CO E mixing ratio of  1 32768 gg E , which is higher than the equilibrium value.

Table 1
List of all Simulations The mixed layer depth can affect convection strength in local-scale cloud-resolving models and may modify the equilibrium climate state (Hohenegger & Stevens, 2016). We use a standard mixed layer depth of 50 m for two reasons. First, it is a reasonable approximation of the typically mixed layer depth in the ocean. Second, a mixed layer depth of 50 m results in a convergence timescale of 5 yr, which we think is fast enough. We examine how different mixed layer depths, including 0.5, 5, 50 and 100 m affect InvCM. The simulations with mixed layer depths of 5, 50, and 100 m take about 0.5, 5, and 10 yr to equilibrate. The fact that the convergence timescale is proportional to the mixed layer depth is consistent with the local SST adjustment being the convergence bottleneck. Interestingly, the simulation using a mixed layer depth of 0.5 m equilibrates at much higher 2 CO E mixing ratio (Figure 1a2). This is associated with spontaneous symmetry breaking, with the highest temperature and precipitation rates occurring off the equator (Figure 2). While this could be an interesting topic of future research, in this article, we focus on simulations with 50 m mixed layer depth.
CO E mixing ratio and the global-mean SST can be linked through ECS for small enough changes. The difference between corresponding SST and the fixed global mean SST (shown in Figure 3c) converges to zero within 5 years. In contrast, normal slab ocean simulations for 310, 320, and 330 K take  45 , 70 , and 65 years, respectively, to equilibrate global-mean SST to within 1 K.
A previous study using a zero-dimensional climate model demonstrates that the e-folding timescale for the surface temperature is proportional to the product of the heat capacity of the system and the temperature timescale is 240 days and the mixed layer depth can be viewed as infinity. Three SST patterns are tested here. The first one is the equilibrium SST pattern in the InvCM-Ensemble group, where the global-mean SST is fixed to 320 K. Then we perturb this SST pattern with to obtain the SST patterns with a lower or higher equator-to-pole SST difference.
change per unit radiative forcing (Cronin & Emanuel, 2013). This timescale is much longer than the SST gradient response timescale because global radiative feedback is the only mechanism that equilibrates the global-mean SST while both atmospheric heat transport and local radiative feedback contribute to adjusting the SST gradient. Since the SST gradient convergence rate is the bottleneck of InvCM (see our Section 2.3), the ratio between the convergence rates of the SST gradient and the global-mean SST is the degree of speedup we gain from the InvCM method. The Cronin and Emanuel (2013) model also explain why the convergence timescale increases with the global-mean SST. As the atmosphere warms, the water vapor content increases nearly exponentially, resulting in the buildup of latent heat in the atmosphere. For hot climates, the atmospheric heat capacity is comparable to, or even greater, than that of a slab ocean with a depth of 50 m. Additionally, the temperature change per unit radiative forcing is large for hot climates (Bloch-Johnson et al., 2015;Bloch-Johnson, Rugenstein, Stolpe, et al., 2020). The long equilibrium timescales we observe in our normal slab ocean simulations could be a result of a combination of the high heat capacity and the large temperature change per unit forcing. Notably, InvCM does not suffer from this phenomenon. In fact, InvCM converges even faster for extremely hot climates (cases of 330 and 340 K, see Figure 3c), which might be due to the fact that increased atmospheric latent heat transport leads to quicker evolution of the local SST.
After equilibrium has been reached for InvCM, the 2 CO E mixing ratio continues to fluctuate, and the amplitude of this fluctuation converted to the corresponding SST depends on the global-mean SST and is largest for the cases of 310 and 320 K (Figures 3a and 3b). The largest amplitude is several Kelvin. This is larger than that of a normal slab ocean simulation in the same conditions and may result in considerable random errors for InvCM. We employ an ensemble of simulations to evaluate such random errors (see the InvCM-Ensemble group in Table 1 for the configuration). Each simulation is run for 10 yr and we take the geometric mean 2 CO E mixing ratio of the last 4 yr as the equilibrium value. Figure 4b shows that the error quantified as the corresponding SST is about 1 K. Only 1 out of 10 members in the ensemble has an error larger than 1 K. For the InvCM-Ensemble group, the global-mean SST is fixed to 320 K, where ECS is highest. A higher ECS is most likely to give bigger errors, so 1K should be an upper bound of random errors for all SST values.
One way to determine the accuracy of InvCM is to compare the equilibrium climate states it produces to those produced by normal slab-ocean simulations. For each simulation in the InvCM-ClimateSensitivity group, we branch a normal slab-ocean run from it. For all branch runs, the deviation of SST is smaller than 1 K ( Figure 5). Some extremely hot simulations exhibit strong variation in SST but show no trend in global-mean SST. How close a state is to equilibrium can also be measured by the TOA energy balance. All branch runs have variations in TOA energy balance less than 1 W 2 m E . This also demonstrates an alternative mixing ratios for several simulations in the InvCM-ClimateSensitivity group. To compare the speed of convergence, the evolution of (d) SST and (e) its difference with the equilibrium SST of simulations in the NormalSlabOcean group are plotted on the same time axis. It takes  5 years for InvCM to equilibrate while the normal slab ocean simulations require 50-70 yr to equilibrate (defined as when the difference between the equilibrium and the current SST is less than 1 K). use of InvCM: it could be used to find states roughly in equilibrium quickly before switching back to a normal slab-ocean simulation. In this way, it would be possible to benefit from the speedup of InvCM and not have to worry about possible systematic bias or a potentially unrealistic configuration.

ECS Bump at High-Temperature Resolution
We investigate the evolution of the climate from 290 to 340 K at an interval of 1 K using 51 simulations (see Table 1, the InvCM-ClimateSensitivity group). First, we focus on the relationship between the equilibrium 2 CO E and SST. Then, we look at the evolution of the distributions of SST, OLR, and ASR. mixing ratio for all 10 members in the InvCM-Ensemble group described in Table 1  mixing ratio is mapped to that in the surface temperature with a scale factor of 20 K, which is the climate sensitivity at 320 K (see Section 3.2 and Figure 6 for climate sensitivity as a function of the global-mean SST).  Table 1, the InvCM-ClimateSensitivity group). We apply a 5 yr moving average to smooth the time-series of the energy imbalance at TOA. can be obtained as a function of the global mean SST (see Figure 6a). We use a 10° polynomial regression to smooth the 2 CO E as a function of surface temperature in the following form: All of the data points lie close to the polynomial fit ( Figure 6a). There is no trend or obvious structure in the residuals and the absolute values of the residuals (no more than  0.05 E doublings of 2 CO E ) are small compared to the overall change in the 2 CO E mixing ratio among the InvCM-ClimateSensitivity group. Moreover, the amplitude of the residuals is of similar magnitude to the random errors in our InvCM-Ensemble simulations (Figure 4b), so these residuals can be accounted for as random errors due to natural variation. We define an infinitesimal climate sensitivity in the following way: where eq E G is the equilibrium 2 CO E mixing ratio as a function of the global-mean SST obtained through a polynomial fit mentioned above. The advantage of this definition is that it defines ECS as a local derivative as opposed to a finite difference, but it still reduces to the usual definition when ECS is constant. In particular: Figure 6b depicts climate sensitivity as a function of SST. The climate sensitivity maximizes near 315 K and the maximum is more than 20 K per doubling of 2 CO E .
The spatial response to the increase of the global mean SST can be captured by InvCM. When the climate warms from 290 to 340 K, the equator-to-pole SST difference decreases significantly (See Figures 7a and 7g). For extremely hot climates (340 K), the equator-to-pole temperature difference is as low as 10K  . Additionally, the SST between  30 E S and 30 N is nearly uniform in this case (Figure 7h). Interestingly, the SST pattern exhibits a small amount of north-south asymmetry for extremely hot climates (Figure 7a). The SST at a latitude of  30 E is always close to the global-mean SST, being only a few degrees higher for a wide range of global mean SST of 50 K. This may not be a coincidence, since half of the area of the planet is equatorward of 30° and half is poleward. Note that the SST can be lower than the freezing point of sea water because we have turned off sea ice for all simulations here. The evolution of the OLR and ASR patterns is more complicated. The OLR at the equator increases dramatically until the global-mean SST reaches 325 K while the subtropical OLR (between  15 E and  30 E ) decreases (Figures 7b and 7c). This trend can be explained by a weakening Hadley circulation. Weaker upwelling at the equator reduces the cloud coverage there and thus results in both higher OLR and ASR. Meanwhile, weaker descent in the subtropics and a warming climate leads to more water vapor and reduces the OLR there (Figures 7d-7f). The positive radiative-convective feedback that Wolf and Toon (2015) described might be the cause of the dramatic decrease in high cloud coverage as SST goes from 310 to 330 K. The global-mean OLR and ASR as a function of SST show a sawtoothed shape with a maximum around 320 K (see Figure 7i). This is similar to the planetary albedo trend reported by Wolf and Toon (2015). There is likely also a clear-sky contribution, as argued by Seeley and Jeevanjee (2021).

Hysteresis and the Onset Runaway Greenhouse
Here, we use a zero-dimensional energy balance model (EBM) to demonstrate the way InvCM could be applied in future GCM studies to fully explore a global climate hysteresis diagram and to investigate climate states very near to the transition to a runaway greenhouse. In our first iteration of the EBM, we define the TOA net downwelling energy flux, ( , ) s E N T G , as: InvCM method of converging a GCM.
First, we consider a climate that can contain hysteresis. The sea-ice feedback has been known to result in hysteresis in planetary climate for a long time (Abbot et al., 2018;Budyko, 1969;Sellers, 1969;Yang et al., 2017), while several recent studies found hysteresis in warm climates due to a positive cloud feedback (Popp et al., 2016;Schneider et al., 2019). Here, we assume that the planetary albedo is specified by the following equation, possibly as a result of changes in stratiform clouds (Schneider et al., 2019): T is the transition temperature scale for the assumed change in stratiform clouds. The smaller w E T is, the more rapid the cloud transition is. If we set w E T = 15 K, there is no climate hysteresis and the equilibrium climate state does not depend on the variable, we choose to evolve (Figure 8a). When we reduce w E T to 10 K, in contrast, climate hysteresis results (Figure 8b). There is now a range of the forcing parameter where three climate states are possible: two stable and one unstable. Critically, we can equilibrate the EBM for all three climate states when we evolve E G (InvCM, Figure 8b2), but we can only equilibrate the EBM for the stable climate states when we evolve s E T (normal method, Figure 8b1). This motivates us to speculate that InvCM could be used to fully specify a climate hysteresis diagram for a GCM by finding unstable climate states.
Next, we study a climate near the transition to a runaway greenhouse state. The runaway greenhouse is a result of the saturation of the negative longwave feedback at high temperatures due to the strong water-vapor feedback (Ingersoll, 2013). We modify our specification of the TOA energy imbalance in our EBM to represent this behavior as follows:

Link to Conventional Climate Modeling
InvCM is an inverse method. Although it does not reveal how the climate responds to forcing directly, it resolves the feedback of the atmosphere due to the change in the global-mean SST and thereby can be used to indirectly find the value of the corresponding forcing. Additionally, InvCM can more easily capture climate states near bifurcations or instability, as pointed out in Section 3.3.
The specific results of normal GCM simulations are not necessarily ends in themselves, in that we do not expect that humanity will exactly quadruple the 2 CO E concentration relative to the preindustrial, or exactly follow any of the emissions scenarios explored in CMIP. Instead, these simulations help us broadly understand the relationship between forcing and response, providing specific instances of a more general mapping. Inverse simulations do the same thing, while running more quickly. In this sense, inverse simulations are no more sensitive to the SST chosen than normal simulations are to 2 CO E concentration chosen-they will both provide similar information about the connection between forcing and response, and allow for similar interpolations to nearby values.
Here, we have implemented InvCM in a relatively simple model configuration (no seasonal cycles and a slab ocean) that is advantageous for theoretical planetary climate studies with GCMs, but unsuitable for realistic near-future climate forecasts. When studying planetary climate, simple climate configurations are often used (Kaspi & Showman, 2015;Komacek & Abbot, 2019;Wolf, 2017;Wolf et al., 2018Wolf et al., , 2020Yang et al., 2013). Moreover, studies sometimes use slab oceans in GCM simulations to investigate the impact of global warming, focusing particularly on the role of the atmosphere (Lu & Cai, 2010;Park et al., 2018). The InvCM method can accelerate similar simulations using a slab ocean. More importantly, the cost of near-global cloud-resolving simulations using a slab-ocean with interactive SST could be drastically decreased using the InvCM method. This is important because although increasingly powerful supercomputers have made near-global cloud-resolving simulations possible, most such simulations currently still use a prescribed SST pattern (Bretherton & Khairoutdinov, 2015;Narenpitak et al., 2017Narenpitak et al., , 2020Tomita et al., 2005).

Possible Extensions and Future Work
In our simulations, the atmosphere is coupled to a global slab ocean with a uniform depth. However, in principle the method we have developed could potentially be modified to incorporate land or sea ice, or to extend the slab ocean to one with spatially varying depths. In this case, it might be advantageous to fix the energy of the ocean instead of the SST. If this is done, the global mean surface temperature may not be constant.
InvCM might also be adapted for a dynamic ocean by applying multiple fluxes at multiple levels. This would address the large heat capacity of a deep ocean, and could lead to a significant decrease in model convergence time. However, the thermohaline circulation is so slow that it can still take thousands of years for a water parcel to close the global overturning circulation, which could limit the decrease in model convergence time. Future research is necessary to determine which effect is more important.
We do not consider seasonal cycles in our simulations. If either the eccentricity or obliquity is nonzero, the seasonal distribution of the shortwave flux will vary. In this case, for normal simulations, the global mean SST may not converge to a constant value, but might instead oscillate in response to the periodic external forcing after it reaches equilibrium. Since InvCM constrains the global mean SST to be constant, it will not typically create the same climate under seasonally varying insolation as a normal simulation. Moreover, since we evolve the 2 CO E mixing ratio to adjust the TOA energy flux to zero, the 2 CO E mixing ratio may oscillate periodically. To reduce the amplitude of this oscillation, we would need to make 2 CO E time scale long compared to a seasonal cycle. As a result, it is possible that InvCM would lose its speed-up advantage, although it would still enable us explore unstable states of an equilibrium curve. Alternatively, if the seasonal forcing is relatively small, we can first run the InvCM simulation until it equilibrates and then uses this state as the initial condition for a normal simulation. Future research is needed to determine the effectiveness of either approach.

Conclusions
We have applied the InvCM equilibration method to the GCM ExoCAM run with a mixed layer ocean. InvCM fixes the global-mean SST and evolves the 2 CO E mixing ratio and SST pattern until an equilibrium is reached. Our main conclusions are: 1. InvCM produces the same climate in ExoCAM as normal slab-ocean simulations with fixed 2 CO E 2. InvCM converges about 10 times faster than normal slab-ocean simulations 3. The speed of InvCM allows us to investigate the climate sensitivity for doubling 2 CO E for global mean surface temperatures ranging from 290 to 340 K at an interval of 1 K. We reproduce the climate sensitivity bump near 315 K at a much higher temperature resolution. This will allow more detailed study of the cause of the climate sensitivity bump 4. InvCM has the potential to find unstable climate states and climate states that are close to the transition to a runaway greenhouse