Revision of “Dependence of Climate States of Gray Atmosphere on Solar Constant: From the Runaway Greenhouse to the Snowball States” by Ishiwatari et al. (2007)

Ishiwatari et al. (2007) (https://doi.org/10.1029/2006JD007368) investigated multiple equilibrium solutions of a gray atmosphere for various values of solar constant, utilizing two types of models, namely, a one‐dimensional energy balance model and a general circulation model (GCM) with simplified physical processes, both of which permit existence of the runaway greenhouse state. The study was retracted, since there was a bug in their GCM that affected quantitative aspects of the study. Here, we revise the study with re‐performing all of the GCM experiments using an appropriately corrected model. The results of re‐experiments show that the main features of the climate regime diagram drawn in the solar constant ‐ ice line latitude plane are mostly unchanged, except that the ice‐free equilibrium state, which existed in Ishiwatari et al. (2007) (https://doi.org/10.1029/2006JD007368) now disappears. It is confirmed that there are the partially ice‐covered state, the globally ice‐covered state, and the runaway greenhouse state. These three states coexist for intermediate values of solar constant. The existence of the large ice cap instability is also confirmed, although the critical latitude of the partially ice‐covered state extends equatorward. Also same as before, the small ice cap instability does not seem to appear.

experiments of INTH2007 by the use of an appropriately revised model with necessary modifications for model tunings. This is the correction study of INTH2007 with the results of those GCM re-experiments. The details of the bug in the model of INTH2007 and the changes from it are described in Supporting Information. Here, we keep, as the substitute of INTH2007, the whole sequence of its arguments and discussions as they were to the extent possible, since we believe that INTH2007 has a certain historical value as the first attempt of systematic investigation of the climate regime diagram of a simple system permitting the runaway greenhouse state and actually there have been some studies referring to INTH2007. We also believe that description of the climatic characteristics of a simple system composed of a gray atmosphere must provide useful insight in considering possible climatic behaviors of more complicated systems. The characteristics of climate with the setups based on the one-dimensional gray atmosphere model of Nakajima et al. (1992) with the swamp ocean is a necessary step for understanding possible diversity of planetary climates.
The following is mostly a repetition of the INTH2007 introduction that is based on the understandings of planetary climate at that time. The major concern of INTH2007 is solar constant that is one of the most important parameters that determine the climate of terrestrial planets. INTH2007 focused attention on the issue that a state of mild climate like that of the present Earth vanishes when the value of solar constant changes by several tens of percent.
When the value of solar constant is decreased from that of the present Earth, the globally ice-covered state emerges. EBMs show that a slight decrease in solar constant results in ice-covered Earth (Budyko, 1969;Sellers, 1969). With the value of solar constant close to that of the present Earth, Budyko's EBM gives multiple equilibrium solutions: a globally ice-covered equilibrium solution and two types of partially ice-covered solutions. The partially ice-covered solution whose ice line is located at the higher latitudes is stable and is regarded as the climate regime of the present Earth. In contrast, the partially ice-covered solution whose ice line is located at the lower latitudes is unstable (Cahalan & North, 1979;Held & Suarez, 1974). This instability is referred to as the large ice cap instability. Note that the number of equilibrium solutions and their stability depend on the setups of EBM. The EBM of Sellers (1969) has yet another unstable equilibrium solution with a small ice cap. The instability is referred to as the small ice cap instability (North, 1984).
Inspired by the snowball Earth hypothesis proposed by Hoffman et al. (1998), GCM experiments on the ice-covered state had been conducted by the time of INTH2007. For example, Baum and Crowley (2001) performed numerical experiments with a coupled atmosphere-ocean GCM and exemplified the possibility of the appearance of the globally ice-covered state under the climate conditions of the Neoproterozoic age. GCM experiments had also been performed to gain insights on the climate on ancient Mars. Abe et al. (2005) explored possible climate states for a wide range of parameters with a terrestrial GCM, and showed that the globally ice-covered state emerges with decreased solar constant regardless of the obliquity of the planet.
When the value of solar constant is increased, the runaway greenhouse state emerges and the atmosphere cannot reach an equilibrium state until the entire ocean evaporates. One-dimensional radiative-convective equilibrium models show that there exists an upper limit to the value of outgoing longwave radiation (OLR) emitted from the top of the atmosphere on a planet with ocean (Ingersoll, 1969;Komabayashi, 1967;Nakajima et al., 1992). When the incident energy flux exceeds this upper limit, no equilibrium state exists. Indeed, experiments with GCMs employing a value of solar constant larger than this limit had demonstrated that the runaway greenhouse state is characterized by the continuous increases both in atmospheric temperature and water vapor amount (Ishiwatari et al., 1998(Ishiwatari et al., , 2002. By the time of INTH2007, it had also been known that one-dimensional radiative-convective models are able to have multiple equilibrium solutions. In the same manner as EBMs, the inclusion of ice-albedo feedback can produce multiple equilibrium solutions (Li et al., 1997). The inclusion of hydrologic cycle also shows the possibility of multiple equilibrium solutions (Rennó, 1997;Sugiyama et al., 2005). For certain values of incoming solar radiation flux, Rennó (1997) showed that two types of stable equilibrium solutions can be obtained: one being an optically thin equilibrium state and the other being an optically thick equilibrium state. Sugiyama et al. (2005) showed that the two types of solutions obtained by Rennó (1997) exist if the derivative of relative humidity with respect to surface temperature is larger than a certain critical value. INTH2007 notified that the introduction of a hydrologic cycle, in general, may change the dynamical ISHIWATARI ET AL.

10.1029/2019JD031761
branch structure of the solutions, although no statistically equilibrium state which could be regarded as that corresponding to the optically thick equilibrium solutions discussed by Rennó (1997) and Sugiyama et al. (2005) was observed in the results of GCM experiments performed by Ishiwatari et al. (2002). Before INTH2007, there had been no attempt to consider the dependence of climate on solar constant by models that permitted both of the runaway greenhouse state and the ice-covered state or in parameter spaces that could cover both of the two states. In the EBMs of Budyko (1969) and Sellers (1969), for instance, OLR is given as a linear function of surface temperature. Hence, there is no upper limit to OLR, and so the runaway greenhouse state cannot be described.
The aim of INTH2007 was to produce the Budyko-Sellers diagram not only with a GCM but with a system setup permitting the runaway greenhouse state. For this aim, the gray atmosphere of Nakajima et al. (1992) was implemented to an EBM and a GCM to examine the branch structure of climatic solutions ranging from the ice-covered state to the runaway greenhouse state. The use of a GCM allowed us to explore climate solutions with explicit representation of the atmospheric circulation and water vapor transport. By comparing the solutions of the EBM and the GCM, it was possible to recognize whether a branch obtained by the EBM had a corresponding branch in the regime diagram obtained by the GCM, and vice versa. It was anticipated that the EBM permitting the runaway greenhouse state might produce a climate regime diagram different from that of Budyko (1969) or Sellers (1969), and that the GCM might produce a climate regime diagram different from that obtained by the EBM.
Since INTH2007, there have been several studies in which GCMs were utilized to investigate climate regime diagrams. The progress in astronomical observations of the last two decades for searching exoplanets has stimulated researches on climates of possible Earth-like planets. For instance, Rose (2015), by the use of an atmosphere-ocean coupled GCM, extended investigation to an ocean planet, that is an aquaplanet but with an explicit model representation of the ocean circulation. Checlair et al. (2017) and Abbot et al. (2018) investigated climate of tidally locked or tidally influenced planets, and Kilic et al. (2017) considered dependence on obliquity. Among those, the research most relevant to INTH2007 is Rose (2015), which shows that the presence of ocean circulation has a significant effect on the branch structure. There appears the water belt regime that is a new stable branch with the ice line in the low latitudinal regions. However, in the results of Rose (2015), the runaway greenhouse state does not appear. The reasons might be the narrow range of the variation of solar constant in addition to the use of a non-gray radiation scheme. Although INTH2007 is now a study over 10 years ago, it seems that the climate regime diagrams including both of the globally ice covered state and the runway greenhouse state has been presented only by INTH2007.
In the followings, we keep the model setups and the structure of arguments as they were in INTH2007. We do not extend our research here, for instance, to implement an ocean circulation model to compare Rose (2015). Further attempt to introduce an ocean model and compare the regime diagram with that of Rose (2015) will be presented in a separate study (Kawai et al., in preparation). The simple setups of the EBM and the GCM utilized here are described briefly in Section 2. The details of the bug in the GCM of INTH2007 and the changes from it are described in Supporting information. The results obtained by the EBM that is the same as that of INTH2007 are re-presented in Section 3. In Section 4, the GCM results are substituted by those re-calculated values and structures by the use of the corrected GCM. The main features of the climate regime diagram drawn in the solar constant -ice line latitude plane with the GCM are mostly unchanged from INTH2007 except that the ice-free equilibrium state now disappears. In section 5, we summarize the results and add some discussion.

Model
The GCM utilized here is DCPAM5 of GFD DENNOU CLUB (http://www.gfd-dennou.org/library/dcpam/ index.htm.en). It has been developed, for the purpose of investigating possible variety of general circulations of atmospheres on terrestrial planets in general (e.g., Noda et al., 2017), by reconstructing AGCM5.3 that was used by Ishiwatari et al. (1998Ishiwatari et al. ( , 2002 and INTH2007. We keep the system setups of INTH2007 to retain present work just as its revision. The framework we adopted is, except for the incorporation of three-dimensional atmospheric motion, basically equivalent to the one-dimensional radiative-convective equilibrium model of Nakajima et al. (1992) in which the existence of radiation limits is clarified and the runaway greenhouse state is described in terms of radiation limits. The adopted configuration of the GCM is the most basic one; we have chosen to retain the simplicity of the model for making an straightforward extension of Nakajima et al. (1992) to a three-dimensional model.
The basic features of the GCM configuration are the same as those of INTH2007 and are briefly summarized as follows. The details of the bug in the model of INTH2007 and the changes from it are described in Supporting Information. The atmosphere consists of a condensable component (water vapor) and a noncondensable component (dry air). Only water vapor absorbs and emits longwave radiation. The absorption coefficient of water vapor is constant and independent of wavelength, 0.01 m 2 kg −1 . Dry air is assumed to be transparent. Cumulus convection is parameterized by a convective adjustment scheme (Manabe et al., 1965). Condensed water is removed from the system immediately without cloud formation. The radiative effects of clouds and the scattering of radiation are excluded. Vertical diffusion is represented by the Mellor-Yamada Level-2 scheme (Mellor & Yamada, 1974). The entire surface is assumed to be a "swamp ocean"; heat capacity is zero and wetness is unity (Delworth & Manabe, 1988). There is no horizontal heat transport by the ocean. Following Budyko (1969), the surface albedo of regions whose surface temperature is below 263 K is 0.5. Otherwise, the surface albedo is zero. During the re-experiment, we found that this stepwise dependence of surface albedo on temperature adversely affects the uniqueness of statistically equilibrium solution. Uniquness of the ice line latitude can be lost since the original formulation of the value of albedo at a grid point is determined only by the surface temperature at that grid point. However, in most part, we present the results with the original stepwise implementation of surface albedo to be adequate for our purpose to substitute the re-calculated results for those of INTH2007. Still, to complete the regime diagram, we present the results with a smoother representation of the surface albedo temperature dependence, where we evaluate the value of surface albedo at each grid point around the ice line latitude considering the proportions of ice-covered and ice-free areas of the "finite surface" represented by that grid point by interpolating the values of surface temperature at the surrounding grid points. Neither snow accumulation, nor ice mass budget, nor sea ice migration is considered. Latent heat of fusion is also excluded from the calculation. This simple model configuration permits the direct comparison of our GCM results with that of EBMs. The dynamical part of the GCM is represented by the pseudo spectral method in the horizontal direction, and by the finite difference method with σ coordinate in the vertical direction. We keep the same spatial resolutions for the re-experiments as those for INTH2007 just to revise the figures and the statistics there. The spectral truncation is the triangular truncation at wavenumber 21 (T21). The number of vertical levels is 16 for the cases where the value of solar constant, S, is smaller than 1450 Wm −2 and 32 for the cases where S is greater than 1450 Wm −2 . The number of vertical levels is increased for the cases with the larger solar constant values, because the tropopause extends to the higher altitudes.
As for the EBM, we did not perform any recalculations. We retain here all of the EBM materials in INTH2007. The EBM utilized there was the same as that of North (1975), except for the radiation scheme. The radiation scheme was the same as that used in our GCM described above, namely, the gray atmosphere of Nakajima et al. (1992). By adopting this radiation scheme, the value of OLR was determined to be a function of the surface temperature on the basis of a one-dimensional radiative-convective equilibrium solution. The relative humidity was given as a fixed value of 60%, which was obtained as a typical value in the GCM results of Ishiwatari et al. (2002). The number of latitudinal grid points for discretization was 181. The meridional heat transport was calculated by the diffusion scheme of North (1975). The value of the heat transport coefficient was 0.5 Wm −2 K −1 . This value was determined in INTH2007 such that the EBM most closely reproduced the meridional temperature difference of the GCM in the ice-free case (fixed zero surface albedo case) with S = 1380 Wm −2 . We assume that the characteristics of the EBM results do not change greatly even if we adapt the values of the EBM parameters tuned to our new GCM results, since the difference between the GCM results of INTH2007 and present one is not very large around S = 1380 Wm −2 and resultant changes of the EBM parameters shall not be very large. In addition to searching for the steady solutions for the EBM, time developments of the solutions were also calculated to investigate the existence of the runaway greenhouse state. In time-dependent problems, the value of heat capacity for a unit area was set to be 1 JK −1 m −2 . The value of heat capacity can be chosen arbitrary, so long as it is not zero, since the value does not affect the results of examination on the existence of the runaway greenhouse state. Zero heat capacity cannot be adopted in the time-dependent EBM, since the tendency term of surface temperature vanishes and surface temperature becomes indefinite.
The value of solar constant S is varied in the range between 1000 Wm −2 and 1710 Wm −2 for the GCM experiments, and in the range between 1200 Wm −2 and 2000 Wm −2 for the EBM experiments. The meridional distribution of incoming solar flux is given as the annual and daily averages of that evaluated using the present terrestrial obliquity but with zero eccentricity: It is symmetric about the equator and varies neither diurnally nor seasonally. As for the GCM runs, the initial condition, if not otherwise mentioned, is the isothermal (280 K) state with constant specific humidity of 10 −3 . In addition, we adopt several other initial conditions, which are rather arbitrarily chosen from calculated results, to pursue the ends of the branches of climate state. The equilibrium solutions of the EBM were obtained by the use of the Powell hybrid method in Minpack which was available at http://www.netlib.org. The runaway solutions of the EBM were obtained by solving the time-dependent problem using an isothermal state having temperature of 330 K or 600 K as the initial conditions. These temperature values were also arbitrarily chosen.

Results of EBM Experiments
The description of this section about the results with the EBM is the same as that of INTH2007. There is no alternation from INTH2007. However, we should remark one thing that we noticed during the GCM re-experiment; the equilibrium solutions which we numerically found for each value of S may not be unique because of the stepwise dependence of surface albedo on temperature and the spatial discretization. The possible lack of uniqueness appears basically in the same manner as described in Held and Suarez (1974) for the EBMs of Budyko (1969) type. Actually, the revised GCM results suffer from this lack of uniqueness as is described in the next section. As for the EBM, we did not pursue this issue since we believe the number of latitudinal grid points for discretization, 181, was fairly large, and the multiplicity of solutions caused by stepwise representation of surface albedo, if any, would affect the structure of the climate regime diagram not greatly. Now, Figure 1 shows the relationship between solar constant and ice line latitude of solutions obtained by the EBM. Ice line latitude at 90° corresponds to an ice-free solution, and 0° to a globally ice-covered solution. Crosses at 90° indicate non-equilibrium solutions where the system is in the runaway greenhouse state.
The gross feature of ice-covered solutions shown in Figure 1 is similar to the results of Budyko (1969) or Sellers (1969). In our EBM, globally ice-covered solutions, indicated as branch α, are found for S ≤ 1903 Wm −2 (point A). One or multiple partially ice-covered solutions are found for 1306 Wm −2 ≤ S ≤ 1903 Wm −2 , indicated by branches β, γ, and δ and critical points A, B, C, and D in Figure 1. Branch β intercepts globally ice-covered solutions (branch α) at point A (S = 1903 Wm −2 ), and branch δ is connected with ice-free equilibrium solutions (branch ɛ) at point D (S = 1436 Wm −2 ). In contrast to the results of Budyko (1969) or Sellers (1969), the branch connections in partially ice-covered solutions and ice-free solutions in the present study display complicated features, the details of which will be explained below.
The branch structure of ice-free equilibrium solutions in Figure 1 displays different features from the results of Budyko (1969) or Sellers (1969). Since the runaway greenhouse state is allowed to emerge in our model, an upper limit to S must exist in order for the ice-free equilibrium solutions to exist. Point E (S = 1441.2 Wm −2 ) is the critical point for the ice-free equilibrium state; no ice-free equilibrium solution can be found in our EBM at any larger values of S. The detail of the branch structure of ice-free equilibrium solution which cannot be recognized in Figure 1 will be shown later, in Figure 2b. For the values of solar constant smaller than that of critical point E, there are two branches of ice-free equilibrium solutions; the first, branch ɛ is connected to partially ice-covered solutions (branch δ) at point D, and the other, branch ζ is connected to the runaway greenhouse state at point F (S = 1310 Wm −2 ). This is the lower limit of S for which the runaway greenhouse state can be found, and it is also the lower limit of S for which an ice-free equilibrium solution can exist. However, the ice-free equilibrium solutions on branch ζ are unstable, as will be explained later in this section.
The partially ice-covered solutions belong to two kinds of branches; one is branch γ in which the ice-covered area decreases with increasing S, and the other consists of two branches β and δ in which the ice-covered area increases with increasing S. As for the stability of the equilibrium solutions, the argument presented by Cahalan and North (1979) should also be applicable for our case since the difference between our radiative term and theirs is not expected to cause a significant difference in the radiative heating property of the partially ice-covered state. The solutions on branch γ are considered to be stable, while the solutions on branches β and δ are unstable. Change of stability occurs at the critical points indicated as point B (S = 1306 Wm −2 ) and point C (S = 1447.5 Wm −2 ). Branch δ has ice line latitudes higher than 69°, which is reached at point C, and branch β has ice line latitudes lower than 28°, which is reached at point B. Following North (1984), the instability corresponding to branch δ whose ice line is located at higher latitudes is referred to as the small ice cap instability, and the instability corresponding to branch β whose ice line is located at lower latitudes is referred to as the large ice cap instability.
In order to clarify the branch configuration of the partially ice-covered state and the ice-free state, the relationship between solar constant S and global mean surface temperature T s of the equilibrium solutions is plotted in Figure 2. Note that the runaway greenhouse state cannot be plotted in these panels, since global mean surface temperature is an increasing function of time for solutions in the runaway greenhouse state. The stable branch γ of the partially ice-covered state corresponds to the branch where global mean surface temperature increases with solar constant in Figure 2a that extends from point B (T s  258K) to point C (T s  298K). Branch β of the large ice cap instability corresponds to the branch extending from point B (T s  258K) to point A (T s  254K). On this branch, when the ice line latitude is far from the equator, global mean surface temperature decreases with the increase in solar constant within the range from point B to S = 1450 Wm −2 . When the ice line latitude becomes closer to the equator, global mean surface temperature increases slightly with an increase in solar constant within the range from S = 1450 Wm −2 to point A. Nevertheless, the equilibrium solutions on this whole branch are locally unstable, as stated earlier.
The configuration of the branch of small ice cap instability is even more complicated. The transition area from the branch of small ice cap instability to the branch of ice-free solutions is enlarged in Figure 2b.  Branch δ of small ice cap instability corresponds to the line segment that extends from point C to point D (T s  297K). It connects to branch γ of the stable partially ice-covered solutions at point C, and to branch ɛ of the ice-free solutions at point D. On branch δ, global mean surface temperature increases with increase in solar constant from point D, reaches a maximum of T s  302K at S = 1442 Wm −2 , and then decreases with solar constant toward point C.
As mentioned earlier, the ice-free state in our EBM consists of two branches, branch ɛ and branch ζ shown in Figure 2b. Point E, The reason for the existence of two branches of the ice-free state in our EBM is best understood from the results of the one-dimensional radiative-convective equilibrium model. It is known that, in a one-dimensional equilibrium model that includes the same radiative process as that of our EBM, there exists two kinds of equilibrium solutions for a given value in a certain range of solar constant (Nakajima et al., 1992). One is the branch of solutions where an increase in solar constant is balanced by an increase in OLR, caused by an increase in surface temperature. The other is the branch of solutions where an increase in solar constant is balanced by an increase in OLR, caused by a decrease in atmospheric opacity due to the reduction in atmospheric water vapor content forced by a decrease in surface temperature. Branch ɛ of ice-free solutions in Figure 2b corresponds to the former branch of the one-dimensional radiative-convective equilibrium solutions, while branch ζ corresponds to the latter. S = 1310 Wm −2 at point F, which is the lower limit of S for the emergence of the runaway greenhouse state, corresponds to the limiting value of OLR for cases with high surface temperature. On the basis of the stability argument on the one-dimensional radiative-convective equilibrium solutions, ice-free equilibrium solutions on branch ɛ are expected to be stable, while those on branch ζ are unstable. Our time-dependent calculations confirm that the ice-free equilibrium solutions on branch ζ are unstable. As is shown in Figure 1, for a given value of incoming solar flux, all equilibrium solutions on branch ζ coexist with a runaway greenhouse solution.

Results of GCM Experiments
We now substitute the re-calculated results by our revised GCM for the results of INTH2007. The major change we found from INTH2007 is the disappearance of the statistically equilibrium ice-free state; there appear runaway greenhouse solutions instead. As for the equatorward edge of the branch of the partially ice-covered state, we had to improve the model implementation of the ice albedo dependence on temperature to obtain a unique solution for the value of solar constant around there.  Figure 1, except for large open circles which indicate the results with the smoother representation of ice albedo temperature dependence. Among them, the upper left one represents the small S end of the runaway greenhouse state, whereas the lower right one, which is superposed by a diamond mark, represents the large S end of globally ice-covered solutions. The remaining ones represent partially ice-covered equilibrium solutions. F represents the results for runs in which the initial condition is the ice-covered solution calculated with S = 1000 Wm −2 , R represents the results for those in which the initial condition is the runaway greenhouse solution obtained with S = 1600 Wm −2 , squares with P represent the results for those in which the initial condition is the partially ice-covered solution calculated with S = 1225 Wm −2 , and asterisks with P represent the results obtained by decreasing S gradually starting from the partially ice-covered solution obtained using S = 1200 Wm −2 . Points with no labels represent the results whose initial conditions are the isothermal state (280 K). crosses. The relationship between solar constant and global mean surface temperature of the statistically equilibrium solutions is shown in Figure 3b. In these figures, results for runs started from various initial conditions are all plotted. Large open circles are the results obtained with the smoother implementation of surface albedo temperature dependence, and will be discussed later.
With the stepwise dependence of surface albedo on temperature (original implementation of surface albedo in INTH2007), partially ice-covered statistically equilibrium solutions are now found for 1150 Wm −2 ≤ S ≤ 1503 Wm −2 . Of these solutions, those denoted by solid circles for 1200 Wm −2 ≤ S ≤ 1502 Wm −2 are obtained from runs with the initial conditions of the isothermal atmosphere of 280K. For S = 1503 Wm −2 , the runaway greenhouse solution is obtained from the isothermal atmosphere of 280K, but, when starting from the partially ice-covered solution for S = 1502 Wm −2 , the partially ice-covered state is maintained (the gray dot labeled with "P"). Further increase of S to 1505 Wm −2 , however, results in the transition to runaway greenhouse state. The solutions for S = 1150 Wm −2 and 1175 Wm −2 denoted by asterisk labeled with "P" are calculated starting from the partially ice-covered solution for S = 1200 Wm −2 by gradually decreasing the value of solar constant. The solution for 1180 Wm −2 denoted by triangle with "R" is calculated starting from the runaway greenhouse state at 6000 days obtained with S = 1600 Wm −2 . A puzzling feature in Figure 3 is that there is another solution for S = 1200 Wm −2 denoted by square labeled with "P," which is obtained from the partially ice-covered solution for S = 1225 Wm −2 . However, we consider the appearance of multiple solutions at S = 1200 Wm −2 does not indicate a new branch of equilibrium solution; it is considered to be caused by stepwise representation of surface albedo. We will discuss and fix the issue shortly below.
Globally ice-covered statistically equilibrium solutions are obtained for S ≤ 1690 Wm −2 . Those for S ≤ 1150 Wm −2 are obtained from the initial condition of the isothermal atmosphere of 280 K. Those labeled by F are obtained from a globally ice-covered solution calculated with S = 1000 Wm −2 . From this initial condition, the globally ice-covered state is sustained up to S = 1690 Wm −2 . At S = 1700 Wm −2 , all of the ice disappears and the runaway greenhouse state appears. The runaway greenhouse state exists for S ≥ 1200 Wm −2 . The solution for S = 1180 Wm −2 starting from a runaway greenhouse solution, as mentioned above, cools down to a partially ice-covered solution. Note also that all of the ice free solutions are in the runaway greenhouse state; a statistically equilibrium ice-free solution obtained for S = 1570 Wm −2 in INTH2007 disappeared.
As mentioned earlier, the multiplicity of solution of the partially ice-covered state with S ≤ 1200 Wm −2 (Figure 3) is considered to be caused by the stepwise representation of surface albedo temperature dependence. Basically the behavior is the same as the loss of uniqueness of the equilibrium solutions for the EBMs of the Budyko (1969) type where the stepwise temperature dependence of albedo is adopted and temperature discontinuity in the solution is permitted at the ice line latitude (Held & Suarez, 1974). Unlike that for the EBMs of Budyko (1969) type, the basic formulation for the GCM assumes temperature to be continuous. However, the discretized model permits an ambiguity of the ice line latitude within the grid interval. We have to use a smoother representation of ice-albedo feedback to fix this issue.
In the smoother implementation of surface albedo, we evaluate the value at each grid point around the ice line latitude considering the proportions of ice-covered and ice-free areas in the "finite surface" represented by that grid point interpolating the values of surface temperature at the surrounding grid points. The re-calculated results with the smoother implementation of surface albedo are plotted as those indicated by large open circles in Figure 3. Now the statistically equilibrium solutions of the partially ice-covered state are uniquely obtained. The equatorward edge of the partially ice-covered state is obtained at around S = 1240 W/m 2 with 20° for the ice line latitude. We confirmed that, this solution jump into the globally ice-covered state when we decrease the value of S to 1230 W/m 2 . The high S edge of the globally ice-covered solution does not change; it is at S = 1690 W/m 2 . If we further increase the value of S to 1700 W/m 2 , then the whole of the ice disappears and the system goes into the runaway greenhouse state. The low S edge of the runway greenhouse state changes a little; the system still is in the runway greenhouse state at S = 1210 W/m 2 , but, when we decrease S to S = 1180 W/m 2 , with the value of which the system with the stepwise representation of surface albedo goes into the partially ice-covered state, it goes into the globally ice-covered state in accordance to the shrinkage of the partially ice-covered state resulting from the change of surface albedo implementation. There is a possibility that the low S edge of the runaway greenhouse state is at slightly larger value of S, but we could not pursue within our computational resource available, since the transition, if any, seems to take very long time.
For S = 1509 W/m 2 , we obtain a partially ice-covered state whose time averaged ice line latitude is 89.4°, where ice disappears frequently but returns repeatedly. The proximity of the ice line to the pole realized there motivated us to explore the existence of ice free statistically equilibrium solutions. We performed a series of long-term GCM integrations for S ranging from 1505 W/m 2 to 1512 W/m 2 without ice-albedo feedback, holding albedo to be that of ice-free surface irrespective of surface temperature. And then, we picked up the solutions that did not go into runaway state, and used them as initial conditions for further integrations but with ice-albedo feedback switched-on. If ice free statistically equilibrium solutions were to materialize, we could expect that the solution should remain to be ice free. The results show that, for all of the ice free initial conditions, ice-covered regions came to develop at least intermittently. Though the longer lasting ice-free period is observed for the larger value of S, there is no completely ice-free solution. The non-existence of ice free statistically equilibrium solution is also confirmed by the temporal evolution of the solution for S = 1510 W/m 2 starting from the partially ice-covered state for S = 1509 W/m 2 (not shown), in which the system evolves to the runway greenhouse state with intermittent appearance of ice in polar region.
The partially ice-covered state extends to the smaller solar constant value and the corresponding ice line latitude reaches lower latitudes compared to those of the stable equilibrium solution for the EBM. The reason why the GCM can maintain the partially ice-covered state at smaller values of solar constant compared to that of the EBM is the existence of the Hadley circulation. In the region of the Hadley circulation, very efficient heat exchange takes place in the meridional direction, and the temperature distribution in the troposphere tends to be latitudinally uniform. In contrast, the heat exchange between the Hadley circulation cells and regions at higher latitudes is comparatively inefficient, and so baroclinically unstable zones with latitudinal temperature gradient are formed (Satoh, 1994). As the value of solar constant decreases, the ice line latitude eventually falls within the Hadley circulation. Then, the equatorward migration of ice line latitude slows down, since the solar flux supplied to the entire tropical region efficiently heats the area of the ice line latitudes because of the efficient thermal mixing in the Hadley circulation. In our EBM, unlike that of Lindzen and Farrell (1977), the latitudinal heat transport is modeled simply as a function proportional to the temperature gradient; thus, the dependence of the efficiency of meridional heat transport on the dynamical difference of the circulation structure is not represented. When the ice line latitude reaches the lower latitudes, heat transport from the equatorial area decreases in the EBM because of the decrease in latitudinal temperature gradient.
INTH2007 mentioned that the structure of the Hadley circulation of the partially ice-covered state with a large ice cap was different from the others; the locations of the precipitation peaks were not at the equator but at the ice line latitudes. Our GCM re-calculations with the stepwise representation of surface albedo also show the similar characteristics but for the solutions with the ice line latitude around 10°, namely, those for 1150 ≤ S < 1200 Wm −2 . They constitute the equatorward end of the partially ice covered state, but they disappear when the smoother implementation of the surface albedo temperature dependence is introduced. Remarkably, the solution at the equatorward edge of the partially ice covered state at around 20° but for the smoother implementation of surface albedo does not show such peculiar characteristics of the Hadley circulation. Nevertheless, we retain the presentation of the typical structure of those solutions, since the atmospheric circulation structure is certainly realized when the ice line is forced to be placed at the latitude, and it has some interesting features. Figures 4 and 5 exemplify the temporal and zonal mean meridional structures of statistically equilibrium solutions of the partially ice-covered state for the case with S = 1380 Wm −2 where the ice line latitude is located at 45°, and for the case with S = 1150 Wm −2 where the ice line latitude is located at around 10°, respectively. These are the results with the same stepwise albedo implementation as INTH2007. The locations of the precipitation peaks for the case with S = 1150 Wm −2 are not at the equator but at the ice line latitudes (Figure 5a), while, for the case with S = 1380 Wm −2 (Figure 4a), the precipitation peak is located near the equator. Corresponding to the distribution of precipitation, the direction of the Hadley circulation in the region of 0.6 < σ < 0.9 for S = 1150 Wm −2 is opposite to that for S = 1380 Wm −2 (Figure 4b). The upward motion in the case of the S = 1150 Wm −2 appears at the ice line latitudes, while downward motion appears at the equator (Figure 5b). Figure 5c shows meridional latent energy transport. Precipitation at the ice line latitudes for S = 1150 Wm −2 is maintained by the latent heat transport from the equatorial latitudes, as shown by the solid line in Figure 5c. The transport of dry static energy is poleward (the solid line in Figure 5d) similar as the ordinary case exemplified for S = 1380 Wm −2 (Figure 4d).
We reconfirmed that the large ice cap instability exists also in the climate system of our GCM. That is, with the smoother implementation of the surface albedo, we could not find a statistically equilibrium partially ice-covered solution having an ice line latitude lower than 20°. There seems to be a critical point around S = 1240 Wm −2 and the ice line latitude of 20°. When we start from the partially ice covered solution for S = 1240 W/m 2 and decrease the value of S to 1230 W/m 2 , the system cannot stay at the partially ice-covered state, but jumps to the globally ice-covered state. On the other hand, the globally ice-covered solutions labeled by F in Figure 3a can be obtained only from a globally ice-covered solution. From an isothermal state of T = 280 K, the system cannot reach the globally ice-covered solutions for S ≥ 1500 Wm −2 , and instead, ends up in either the partially ice-covered state or the runaway greenhouse state. These results suggest that an unstable branch is present in the GCM solutions connecting the small S end of the partially ice-covered state branch and the large S end of the globally ice-covered state branch.
As for the small ice cap instability, the transition from the partially ice-covered to the ice-free states appears to be continuous, being similar to that observed in INTH2007, except that the ice-free state is not a statistically equilibrium state but the runaway state in the present re-experiment. In short, the small ice cap instability does not exist. The smooth transition from the partially ice-covered state to the ice-free state corresponds to the increase of temporal variability as shown in Figure 6. As the value of solar constant increases and the ice line latitude recedes poleward, the fluctuation amplitude of the ice line latitude increases ( Figure 6, note these are the results by the stepwise representation of surface albedo implementation). This result seems to be similar to the results of Lee and North (1995) where the continuous transition from partially ice-covered solutions to the ice-free solutions are obtained with a GCM and a noise-forced EBM. In our re-experiment, however, the ice-free solutions are in the runaway greenhouse state. Figure 6 also seems to contain a reason why the branch of the partially ice-covered state is uniquely defined for the range of larger S even when the stepwise representation of surface albedo is used. When the value of S is large, Figure 6 shows the fluctuation amplitude of the instantaneous ice line latitude increases. Resultantly, the temporal mean ice line latitude, which is the latitude indicated in the regime diagram Figure 3, becomes more naturally represented as a latitudinally continuous quantity in spite of discretization, and is allowed to move smoothly with the change of S. When the value of S is small, on the other hand, the fluctuation amplitude of the instantaneous ice line decreases, and once the amplitude becomes smaller than the model resolution, the discontinuous nature of the albedo representation becomes surfacing and the ice line tends to stick to the grid points of the GCM. ISHIWATARI ET AL.   The runaway greenhouse state of the GCM emerges, when the initial condition is the 280 K isothermal state, for S ≥ 1503 Wm −2 and for S ≥ 1510 Wm −2 with the stepwise and the smoother surface albedo implementations, respectively. When a runaway greenhouse solution is used as the initial condition, the range of S which maintains the runaway greenhouse state extends to the smaller S values. These correspond to the solutions marked with label R in Figure 3a, for which the runaway greenhouse solution obtained with S = 1600 Wm −2 is given as the initial condition. The runaway greenhouse state is maintained even when the value of solar constant is reduced to S = 1200 Wm −2 (to S = 1210 Wm −2 with the smoother surface albedo implementation); globally averaged surface temperature continues to increase. However, for the value of solar constant S = 1180 Wm −2 , surface temperature begins to drop; the globally ice-covered state is obtained with the smoother surface albedo implementation. The smallest value of S which can maintain the runaway greenhouse state can be roughly described with the limiting value of OLR which is achieved in the case where the atmospheric water vapor amount is sufficiently increased in the one-dimensional radiative-convective model having a fixed value of tropospheric relative humidity (Nakajima et al., 1992). The asymptotic value of OLR, which is realized in the limit of a water vapor atmosphere, estimated by one-dimensional model is about 300 Wm −2 . This estimation is in fairly good agreement with the lowest limit of the globally mean incoming solar flux that brings about the maintenance of the runaway greenhouse state (S/4).

Summary and Discussions
We have revised INTH2007 reflecting the results of experiments re-performed with the appropriately corrected GCM. The major change we found is the disappearance of the statistically equilibrium ice-free state; there appear runway greenhouse solutions instead. We also found that the stepwise dependence of surface albedo on temperature should be implemented carefully. Other major results of INTH2007 did not change qualitatively. The following "Summary and Discussion" is mostly the same as that of INTH2007 with a few quantitative changes of critical values of solar constant.
The results obtained by our GCM suggest that the large ice cap instability also occurs in a three-dimensional system. We may speculate that there exists an unstable solution with a large ice cap whose ice line latitude is located between 20° and the equator. Such an unstable solution in the GCM, if it exists, will be difficult to find. The structure of the large ice cap instability of the GCM may be altered from that of the EBM because of the existence of Hadley circulation in the GCM. The EBM of Lindzen and Farrell (1977), which takes into account the effect of Hadley circulation, may be useful in examining the contrast between the large ice cap instabilities of the GCM and the EBM in detail. However, we leave this theme for future studies. As for the small ice caps, as in INTH07, our GCM continuously yields statistically equilibrium states up to the value of S for which the ice-line very close to the pole results, which seems to indicate that, in our three-dimensional system, the small ice cap instability does not exist. Only a notable difference from INTH07 is that the branch of the partially ice covered state is directly connected to the runaway state in the present re-experiment instead of the ice-free statistically equilibrium state. Since Lee and North (1995) predicts that the manifestation of small ice cap instability depends on the strength of noise forcing, a GCM experiment including the effect of the response time of ice cover might be an interesting subject. The thermal inertia of ice sheets may affect stability. In our undergoing experiment, which has actually triggered this revision, it is found that existence/non-existence of ice-free equilibrium states depends on specifications of heat capacity of the ocean and heat transport by the ocean. This problem will be discussed in the separate study (Kawai et al., in preparation).
On the results of the EBM, we should remark one thing that we noticed during the re-experiment of the new GCM; the equilibrium solution which we numerically found for each value of S may not be unique because of the stepwise dependence of surface albedo on temperature and the spatial discretization. The possible lack of uniqueness appears basically in the same manner as described in Held and Suarez (1974) for the EBMs of Budyko (1969) type. Actually, the GCM results with the stepwise representation of surface albedo temperature dependence suffer from this lack of uniqueness as is described in Section 4. As for the EBM, we did not pursue this issue since we believe the number of latitudinal grid points for discretization, 181, was fairly large, the multiplicity of solutions caused by stepwise representation of surface albedo, if any, will affect the structure of the climate regime diagram not greatly.
The statistically equilibrium states obtained by our GCM correspond to some of the equilibrium solutions appearing in the vertically one-dimensional radiative-convective models. The corresponding ones are the equilibrium solutions on the branch of increasing OLR with increase in surface temperature obtained by Nakajima et al. (1992), and the optically thin equilibrium solutions of Rennó (1997) and Sugiyama et al. (2005). The GCM results of this study or of Ishiwatari et al. (2002) do not seem to show any correspondences to the optically thick equilibrium solutions of Rennó (1997) and Sugiyama et al. (2005). The difference lies in the dependence of relative humidity on surface temperature. In the optically thick equilibrium solutions of Sugiyama et al. (2005), relative humidity increases with increase in surface temperature, while, in the tropics of our GCM, relative humidity decreases (Ishiwatari et al., 2002). There is a possibility that the moistening efficiency of cumulus might affect the branch structure of solutions. Rennó (1997) mentions that the optically thick equilibrium solutions could not be found when convective adjustment scheme was utilized or when relative humidity was fixed in his one-dimensional radiative convective model. Solutions corresponding to the optically thick equilibrium solutions of Rennó (1997) and Sugiyama et al. (2005) might be obtained in the GCM if we use a cumulus parameterization scheme that has a vigorous moistening effect against the drying in the Hadley circulation region. However, it is unclear at this point whether the moistening effect of cumulus would be so efficient when surface temperature is increased.
Our GCM demonstrates that a runaway greenhouse state or a statistically equilibrium state can be realized for the same value of solar constant in the range of 1240 ≤ S ≤ 1509 Wm −2 with the smoother surface albedo implementation. However, we have not been able to confirm the expectation of Rennó (1997) that, in order to produce a runaway greenhouse regime from an equilibrium solution, the finite amplitude perturbation must be large enough to produce a surface temperature larger than a critical value which corresponds to the unstable solution of Nakajima et al. (1992). The initial perturbation utilized here for activating the runaway state for 1210 ≤ S ≤ 1509 Wm −2 in our GCM is quite huge; it is the runaway solution at S = 1600 Wm −2 . At the moment, we have not been able to determine a minimal magnitude of perturbation necessary to produce a runaway greenhouse solution from an equilibrium solution. In other words, it is required to figure out the GCM equivalent for the unstable equilibrium solution of Nakajima et al. (1992).
The critical values of solar constant at which the number of equilibrium state changes as shown in Figure 3a may vary according to the model configurations, such as the value of albedo, the details of the radiation scheme, and especially, the radiation effect of cloud which is excluded in the present study. However, the importance of our results is that we have exemplified the possibility of the coexistence of the globally ice-covered state, the partially ice-covered state, and the runaway greenhouse state, for a given, common solar constant value. The actual state to be realized at a particular value of solar constant depends on the initial condition and how the value of solar constant is varied. Our result shows that the runaway greenhouse state, which cannot be represented by previous EBMs (e.g., Budyko, 1969), is located not very far from the partially ice-covered or the globally ice-covered state in the climate regime diagram. Although the results of those previous EBMs are usually referred to as the standard climate regime diagram, it should be noted that their diagrams are missing an important regime: the runaway greenhouse state.
In this study, a climate regime diagram for an aquaplanet with gray atmosphere is presented aiming as a straight forward extension of Nakajima et al. (1992). The next step toward a general understanding on the variety of climates of terrestrial planets is to examine the climate regime diagram with implementing the oceanic circulation, like that performed by Rose (2015), but for this gray atmosphere. On this issue, we are now preparing the separate study Kawai et al. (in preparation) which will provide information about the effects of dynamic oceanic on the regime diagram and discussion about the difference from the diagram of Rose (2015).

Data Availability Statement
Model source code and configuration files necessary for numerical experiments are available at https://doi.org/10.5281/zenodo.3514109.