Two‐Dimensional Idealized Hadley Circulation Simulation for Global High Resolution Model Development

Abstract

A new framework for the development of global high resolution models using a vertical‐meridional, two‐dimensional pole‐to‐pole domain is proposed. Compared with the three‐dimensional global model, the two‐dimensional framework can simulate cloud‐scale convection at significantly reduced computational cost, data volume, and turn‐around time, and accelerate the development phase. Although idealized, it allows for simulation of various clouds and scale interactions under a wide range of environmental conditions, even though zonal waves such as mid latitude storms cannot be simulated. Physics schemes can be tested and assessed over a full array of simultaneously occurring cloud regimes. The framework can serve as an intermediate step between a one‐dimensional simulation and a full three‐dimensional global model. Using this framework, we analyze the resolution dependencies of simulated clouds and the Hadley circulation in a range of global high resolution models by performing multiple 1000‐day simulations; for some resolutions tested, for example, 2 km horizontal resolution, three‐dimensional simulation is not feasible at the present time. Analysis shows that both horizontal and vertical resolution determine the properties of the Hadley circulation and clouds, and that the differences among the simulated Hadley circulations are understood as differences in scale interactions influenced by resolution.

Introduction

The field of atmospheric science is heading towards an era of global high resolution modeling, which uses a horizontal grid spacing of O(1–10 km) and will be used to assess climate projections in conjunction with conventional low resolution global circulation models (GCMs) that nowadays typically operate with 1–2° resolution. In a global high resolution model (GHM), also referred to as a “global cloud system resolving model” or “global storm resolving model,” convection is not parameterized but is represented explicitly by a fluid dynamics core and physics parameterizations that represent microphysics, turbulence, and radiation (e.g., Hohenegger et al., 2020; Miyamoto et al., 2013). Recently a GHM intercomparison project with 5 km grid spacing was carried out to explore performance improvement using grid refinement (DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains, DYAMOND; Stevens et al., 2019). In the project, nine participating models showed very realistic clouds compared to satellite observations. Moreover GHM is now starting to be used in weather forecasts and climate studies (Bretherton & Khairoutdinov, 2015; Narenpitak et al., 2017; Wedi et al., 2020).

Improvement in GHM performance with finer grid spacing is mathematically reasonable in a spatially discretized system. However, resolution improvement cannot solve all of the problems simultaneously and automatically since assumptions and parameters are still included in the parameterizations for microphysics, turbulence, and radiation—that is, the performance and the behavior of physics schemes are influenced by the resolution. Since these physics schemes in GHMs are expected to show different manifestations from those used in conventional GCMs, which rely on convective parameterization, (e.g., Dueben et al., 2020), they should be replaced by new schemes, or tuned for a target resolution. Thus, to demonstrate their applicability, performance assessments for new schemes or new parameter sets need to be performed at the resolution to be used in production runs. Unfortunately, performing multiple GHM simulations in the model development phase is computationally too expensive at the present time, and is exacerbated by enormous volumes of data.

Single column models or regional models have often been used in the development phase to avoid large computational cost and to accelerate the parameter adjustment/tuning process. Vertical resolution can be increased easily in a single column model with the caveat that performance related to the horizontal variability cannot be assessed with a single column model. Although a three-dimensional regional model can be run at high resolution and incorporate horizontal variability, the simulated interactions between the large-scale and the small-scale is limited due to the absence of a simulated large-scale circulation. Additionally, the computational cost becomes large when the resolution is extremely high. Although single column models and regional models are good options for the initial stage of development, they do not enable one to evaluate horizontal structures and/or scale interactions at an advanced stage of the development. In this advanced stage, a coarse resolution GCM has frequently been used to evaluate new physics schemes. When the resolution used in the production run is not very different from that used in the developmental test runs, the evaluation is anticipated to work well. Compared with a grid spacing of a GHM production run, however, a grid spacing of 1° or 0.25° is too large and the comparison may not be able to address the model evaluation, simply because a performance problem found in a coarse resolution GCM may be different from that found in a GHM.

We propose a vertical-meridional two-dimensional “global” simulation as a framework to accelerate model physics development and tuning for GHMs. Here the two-dimensional “global domain” is a domain extending from the south-pole to the north-pole and including the entire troposphere. By reducing the zonal dimension from the three-dimensional domain, computational cost and data volume can be reduced dramatically. This allows one to use any desired resolution for storm resolving simulations and to perform numerous test cases without demanding extremely large computational resources, as in the case of a three-dimensional high resolution global simulation. Performance assessments also benefit from fewer degrees of freedom in the two-dimensional domain than the three-dimensional domain; one expects a converged numerical solution and convergence behavior to be more straightforward in a system with fewer degrees of freedom. Another advantage of this framework is that one can study behaviors of the new physics schemes in a balanced state by running simulations for sufficiently long duration.

The usefulness of the two-dimensional framework for convection and cloud system in cloud system resolving model has been demonstrated by previous studies. Grabowski (1998) and Wu et al. (1998) showed that two-dimensional simulations can be used as realizations of deep convection and cloud systems in the Tropics. Grabowski (2001) applied a two-dimensional cloud-resolving model to represent cloud processes not only for the Tropics but also for the Subtropics and the Mid-latitudes.

A caveat of the two-dimensional global simulation framework is a lack of observations for comparison. Thus, an important question to be considered is whether trends and relationships among important variables observed in the two-dimensional framework are in fact similar to these in GHMs, allowing the two-dimensional framework to serve as an intermediate step between single column model and GHM. During the course of the paper, we will provide examples of common behaviors and trends between our two-dimensional simulations and the existing literature for GHM simulations to support the usefulness of the two-dimensional framework. Once consistency is established among two-dimensional global models, and various shared behaviors and trends between two-dimensional global models and GHMs are confirmed, the knowledge acquired will provide guidance as to how one might develop a GHM at one's target resolution.

As shown below in a two-dimensional global simulation, an idealized Hadley circulation, which is a good representation of Earth's circulation, can be reasonably simulated on the vertical-meridional plane. The climate is ultimately a result of an energy balance between solar radiation input and energy output from the Earth, and the Hadley circulation is the primary response of the Tropical atmospheric general circulation to the global temperature distribution (E. K. Schneider, 1987; Oort & Yienger, 1996). In our two-dimensional global simulation, a prescribed sea surface temperature (SST) imposes the meridional distribution of the temperature gradient instead of the incoming solar radiation. A large-scale convective circulation is organized as a response to the meridional distribution of the temperature gradient in each hemisphere, corresponding to the Hadley circulation in the real atmosphere. The idealized Hadley circulation in the vertical-meridional plane can be considered as a steady state of the Hadley circulation in a three dimensional model (E. K. Schneider, 1987). In other words, the vertical-meridional two-dimensional simulation is equivalent to simulating a zonal mean circulation without eddy diffusion.

In this study, we report on resolution sensitivity tests as a useful example of an application of this two-dimensional global simulation framework. Our motivation is to investigate how the large-scale circulations and clouds embedded in the circulation change due to changes in horizontal resolution (2–16 km) and vertical resolution (128 and 256 levels). We show results obtained from 19 1000-day two-dimensional global simulations to investigate a quasi-balanced state. This type of sensitivity test at high resolution and for long duration is almost impossible for three dimensional GHMs considering current computer power and computational cost. Hohenegger et al. (2020) studied the effects of horizontal resolution on a GHM for a horizontal resolution range from 80 to 2.5 km over a simulation duration of only 40 days. Dueben et al. (2020) investigated performance differences by comparing a GHM at 1.45 km horizontal resolution to a conventional GCM at 9 km horizontal resolution in simulations of only 12 or 24 hr. Their simulations are too short to examine the dependence of cloud properties in a balanced state and large-scale circulation on resolution. Prein et al. (2015) demonstrated the advantages of high resolution simulation using regional climate models at a sub-10 km resolution, and compared to conventional GCMs at 1 degree resolution. Dueben et al. (2020) and Prein et al. (2015) include the differences between explicit convection and the parameterized convection in their comparison, which makes it difficult to discern an effect from resolution alone. The current work will explore these issues more systematically, and for longer integration times.

The rest of this study is arranged as follows. A detailed description of the vertical-meridional two-dimensional framework is given in Section 2. The general characteristics of the simulations are described in Section 3. We discuss resolution dependencies for the simulated large-scale circulation and clouds in Section 4 and possible mechanisms for the resolution dependencies in Section 5. In Section 6, a summary of the study is presented.

Numerical Model and Experimental Settings Governing Equations for a Two-Dimensional Atmosphere, and Model Physics

The numerical model is based on the System for Atmospheric Modeling (SAM, version 6.10.10; Khairoutdinov & Randall, 2003). SAM consists of an anelastic governing equation using the finite difference approximation formulated on the Arakawa C-grid with a height coordinate. Velocity components are predicted using the third-order Adams-Bashforth scheme (Durran, 1991) and the second-order center advection scheme. Scalar variables (liquid water static energy, subgrid scale turbulence kinetic energy, and prognostic microphysical variables) are predicted with the forward in time, monotonic fifth-order advection scheme of Yamaguchi et al. (2011). Smolarkiewicz et al. (2001) demonstrated an application of the anelastic governing equation to a global simulation. They concluded that the anelastic global model can simulate atmospheric flows adequately, while minor differences were found from other non-hydrostatic governing equations. Bretherton and Khairoutdinov (2015) also demonstrated a practical application of SAM for global scale simulations.

We modified the momentum equations of SAM to achieve angular momentum conservation on the vertical-meridional, two-dimensional plane following Satoh (1994) using the shallow water approximation and neglecting the cosine Coriolis force and the curvature terms. The governing equations for the continuity equation, liquid/ice static energy, and hydrometeors are not modified. In the modified governing equations, the prognostic variables for momentum are meridional velocity (v), vertical velocity (w), and angular momentum (l). The velocity in the zonal direction (u) is diagnosed through the angular momentum (Equation 3 below) considering conservation of angular momentum. [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where t is time, x_i (i = 1, 2) are meridional and vertical directions (i.e., y and z), and u_i (i = 1, 2) are the wind components along y and z directions (i.e., v and w). ρ is air density and p′ is perturbation pressure relative to the reference state. τ_ij is the sub-grid scale (SGS) stress tensor, which is given by the SGS turbulence scheme. F_l is the SGS angular momentum flux. B is the buoyancy term and δ is the Kronecker delta. Ω, a, and φ are the angular velocity of Earth's rotation, radius of the Earth, and latitude, respectively. The Coriolis factor f is defined by Equation 5.

Angular momentum, defined by Equation 3, is treated as a scalar variable in the model. The meridional velocity includes the Coriolis force in the momentum equation. In this framework, the angular momentum works as a kind of external forcing to the meridional velocity through the diagnosed zonal velocity. The angular momentum has source and sink terms only at the bottom boundary where a surface exchange of the angular momentum is diagnosed through surface wind drag in the zonal direction, defined by Equation 4 based on a bulk formula. C_d is the drag coefficient, V is the magnitude of the surface wind, and U₁₀ is the zonal velocity corresponding to 10 m height. For pressure, advection, and diffusion terms, we modified the original schemes available in SAM by removing the calculations of partial derivatives in the zonal direction, so that these terms are calculated on the vertical-meridional two-dimensional plane. The minimum grid number in the zonal direction is four, as required by SAM's infrastructure. We simply copy one grid point to the other three grid points in the zonal direction to ensure two-dimensionality and eliminate the degree of freedom in the zonal direction.

We use a two moment bulk microphysics scheme, the Predicted Particle Properties (P3; Milbrandt & Morrison, 2013; Morrison & Milbrandt, 2015) for cloud microphysics process; the Rapid Radiative Transfer Model for GCMs (RRTMG; (Iacono et al., 2008; Mlawer et al., 1997)) for shortwave and longwave radiation processes; and Monin-Obukhov similarity (Monin & Obukhov, 1954) for surface fluxes.

The 1.5-order turbulence kinetic energy (TKE) scheme based on (Deardorff, 1980) with explicit diffusion scheme is the default turbulence scheme in SAM, which hereafter we refer to as TKE1.5. Cheng et al. (2010) and Bogenschutz and Krueger (2013) argued that TKE1.5 cannot simulate various scales of convection adequately at the typical grid spacings used for cloud system resolving simulations. For this reason, we employ the Simplified High Order Closure (SHOC; Bogenschutz & Krueger, 2013) turbulence scheme instead of TKE1.5. Our analysis of SHOC, however, reveals that a combination of SHOC and SAM's default explicit diffusion scheme tends to blow up with grid spacings of several km due to SHOC being a non-local mixing scheme. To avoid numerical instability, we modified SHOC so that it behaves like TKE1.5 for small grid spacings, and like the original SHOC for large grid spacings. Details of the modification and the performance checks are presented in Appendices A and B. In addition, we use only the turbulence closure part of SHOC, even though SHOC provides diagnosed SGS cloud fraction and SGS cloud mass, as discussed in Appendix B. Thus, condensation is computed in P3 assuming saturation adjustment.

Experimental Settings

The computational domain is a two-dimensional plane that covers the meridional direction from 90°S to 90°N and the vertical direction from the surface to 35 km height. Horizontal grid spacings are 2, 4, 8, and 16 km and the corresponding numbers of grid points are 10, 240, 5,120, 2,560, and 1,280, respectively. Eight types of vertical grid arrangements are used (Table 1). Hereafter, the notation convention given in the table is used for each vertical grid type. The six vertical grids, L128-DZ35, L128-DZ50, L128-DZ100, L256-DZ30, L256-DZ50, and L256-DZ100, are constructed as follows: (a) the number of levels is either 128 or 256, (b) a constant, minimum grid spacing (Δz) (30, 35, 50, or 100 m) is used up to 2 km height, (c) the grid spacing from 2 km height to 15 km height is the same to the extent possible amongst each vertical grid group (either L128 or L256), and (d) the grid spacing is stretched smoothly from a height of 15–35 km (domain top). L128-DZ50z1 and L128-DZ100z1 are the same as L128-DZ50 and L128-DZ100 but with a lowest layer thickness of 35 m to investigate the impact of the first layer. The number of vertical levels used in the DYAMOND project (Stevens et al., 2019) ranged from 74 to 137, but the model top height was 80 km for 137 vertical levels. The 128 vertical level setup in this study is similar to that in the DYAMOND project and the 256 vertical level grid has approximately twice the largest number of levels in the DYAMOND project.

Table 1 Vertical Grid Arrangement: The Columns for Stretched Δz Show the Mean (Maximum) Value in the Levels

Notation convention	Grid number	Constant Δz bottom to 2 km	Stretched Δz 2 km to 15 km	Stretched Δz 15 km to 25 km
L128-DZ35	128	35 m	259 m (400 m)	747 m (1,175 m)
L128-DZ50	128	50 m	271 m (400 m)	448 m (503 m)
L128-DZ100	128	100 m	308 m (400 m)	400 m (400 m)
L128-DZ50z1	128	50 m (35 m at z = 1)	271 m (400 m)	450 m (509 m)
L128-DZ100z1	128	100 m (35 m at z = 1)	308 m (400 m)	400 m (400 m)
L256-DZ30	256	30 m	116 m (150 m)	181 m (367 m)
L256-DZ50	256	50 m	124 m (150 m)	157 m (190 m)
L256-DZ100	256	100 m	138 m (150 m)	150 m (152 m)

Note. The information above 25 km is omitted from the list, because the levels above 25 km are in a damping layer.

Each simulation case is named by a combination of horizontal grid spacing in a format of DY# and the type of vertical grid; for example, a 4 km horizontal grid spacing with 128 vertical levels with the minimum Δz = 35 m is referred to as DY4-L128-DZ35. All 19 simulation cases performed in this study are listed in Table 2. Hereafter, we use “high resolution” to convey “fine grid spacing” and vice versa. Since the difference in vertical resolution between DY4-L128 and DY4-L256 cases is basically in the free atmosphere (FA), we can evaluate the FA-resolution impact by comparing these two groups. Similarly, the difference in vertical resolution among DZ35, DZ50, and DZ100 cases is basically in the planetary boundary layer (BL), so we can assess the effect of BL-resolution using these three vertical grid arrangements.

Table 2 The List of Simulation Cases: Cross Marks Indicate the Cases Carried Out in This Study

	2 km	4 km	8 km	16 km
L128-DZ35	×	×	×	×
L128-DZ50		×	×	×
L128-DZ100		×	×	×
L128-DZ50z1		×	×	×
L128-DZ100z1		×	×	×
L256-DZ30		×
L256-DZ50		×
L256-DZ100		×

The lateral boundary, that is, the north and south poles, is periodic with buffer regions around both poles in order to reduce instabilities induced near the poles. The buffer is implemented as a Newtonian cooling, and damps meridional and vertical velocities to zero. The width of the buffer region is 64 grids (32 grids for 16 km grid spacing) from each lateral boundary. The damping factor is reduced toward the inner domain. Furthermore, in order to reduce wave reflections from the model top, damping for meridional and vertical velocities is applied to the upper 30% of the vertical levels above 25 km.

The bottom boundary condition is ocean with a fixed SST profile for simplicity, and sea ice is not taken into account. The SST profile is derived from the ECMWF Reanalysis v5 (ERA5; Hersbach et al., 2020) by averaging temporally over 12 years from 2000 to 2011 and zonally from 0° to 360°. Figure 1 shows the SST profile for the Northern Hemisphere. The maximum value of the SST is 301.5 K at 6.5°N. For a simple, idealized simulation, the SST profile from the Northern Hemisphere is copied to the Southern Hemisphere: the SST profile is identical for both hemispheres and there are two maxima in the Tropics, that is, a double inter-tropical convergence zone condition.

Figure 1. Meridional profile of the temporal and zonal mean sea surface temperature derived from ECMWF Reanalysis v5.

Each simulation is started from a steady state: the initial conditions for v and w are zero. The initial condition for angular momentum is calculated for the steady state (i.e., u = 0), and the initial surface pressure is set as 1,013 hPa for the entire domain (U.S. standard atmosphere (NOAA, 1976)). The pressure profile in each column is calculated assuming the hydrostatic balance using the surface pressure and air temperature given by ERA5. The two physical parameters used in our modification are set as in the real atmosphere; Ω = 7.27 × 10⁻⁵ rad s⁻¹ and a = 6.37 × 10⁶ m. The aerosol number concentration in the P3 microphysics scheme is fixed at 100 mg⁻¹. The other parameters are the same as the default setting in SAM. The simulation duration is 1,000 days with a 10 s time step for all cases. The time step is automatically divided into sub-steps when the Courant-Friedrichs-Lewy condition is anticipated. The simulation date is fixed at the autumn equinox day for the duration. The diurnal cycle of solar insolation is included. These configurations for physics schemes are identical for all simulation cases.

While six simulation cases with the vertical grid arrangement of L128-DZ50z1 and L128-DZ100z1 were carried out as shown in Table 2, we did not find a major impact from the 35 m height of the lowest model layer. We briefly summarize the results for L128-DZ50z1 and L128-DZ100z1 cases in Section 5. In Sections 3 and 4, we discuss the remaining 13 cases listed in Table 2.

We also conducted nine simulation cases with the TKE1.5 scheme for the combinations with 4 km, 8 km, and 16 km horizontal grid spacings and with L128-DZ35, L128-DZ50, and L128-DZ100 vertical grid arrangements. Thus, a total of 28 simulations were performed for a range of storm resolving resolutions. The simulation results with the TKE1.5 scheme are shown in Appendix D.

The rundown for numerical cost and data size are as follows: A 1000 days time integration for the DY4-L128-DZ35 case costs 100 hr of CPU time using 256 CPU cores (Intel Skylake) with 8 nodes of the NOAA Research and Development HPC System supercomputer and 760 hr of CPU time using 512 CPU cores (Intel Knights Landing) with 8 nodes of the DOE National Energy Research Scientific Computing Center supercomputer. In the horizontal 4 km grid spacing simulation, approximately 0.6 simulated-years per day is achieved using only eight nodes of a supercomputer. The total data size of the two-dimensional history output for a single case is 113 GB including 23 variables for 2001 snapshots with an output frequency of twice per a day.

General Characteristics of the Simulated Hadley Circulation

Common characteristics of the simulated large-scale circulation and clouds are described in this section. Figure 2 shows the temporal evolution of the global accumulated angular momentum normalized by the initial value for various combinations of horizontal and vertical grid spacings. Since the initial condition is a steady state, the initial value is 1.0 for all cases in Figure 2. The global accumulated angular momentum increases rapidly after the initiation of the simulation, and reaches a quasi-balanced state around the 300th day; fluctuations can be observed thereafter, with a period of about 100 days. The global accumulated angular momentum in the quasi-balanced state is approximately 1.03, which indicates that the mean flow is westerly. We define an analysis period for the quasi-balanced state as the last 700 days (i.e., from the 300th day to the 1000th day) based on this temporal evolution of the global accumulated angular momentum.

Figure 2. Temporal evolution of global accumulated angular momentum normalized by the initial value. The gray shading shows the range between maximum and minimum for all 19 cases.

Figures 3 and 4 show vertical-meridional cross sections averaged over the analysis period and both hemispheres for the DY2-L128-DZ35 case. Since the prescribed SST is symmetric about the equator and the date of the simulations is the autumn equinox, the 700-day averaged circulation is almost identical for both hemispheres. Thus, a hemispheric latitude-height plane is used to show the results in the figures in this study. A two-dimensional mass-flux stream function is calculated by vertically integrating the mean meridional wind (Mitas, 2005; Oort & Yienger, 1996; Waliser et al., 1999), [Image Omitted. See PDF]where p_t and p_s are pressure at the top and the bottom of the model respectively; the overbar means temporal mean about the analysis period and g is gravitational acceleration. We assume that Ψ = 0 and v = 0 at the model top. Positive Ψ corresponds to a clockwise circulation and a negative Ψ to a counter clockwise circulation on the vertical-meridional plane (Figure 3a). Note that the circulation direction of the mass-flux stream function in the southern hemisphere is flipped (i.e., multiplied by − 1) when calculating the hemispheric mean. In Figure 3a, the red colored area, that is, positive Ψ and clockwise circulation, ranging from 10° to 45° corresponds to the Hadley circulation. A counter circulation indicated by blue colors, that is, negative Ψ and counter clockwise circulation, is found in the equatorial area, which implies that the simulated Hadley circulation has a structure like the double inter-tropical convergence zone. Recall that the computational domain spans from 90°S to 90°N and that there are two SST maxima at around 6.5° in both hemispheres where the strongest convection tends to be organized (Figures 4a and 4d). Thus, this counter circulation is not an artificial structure, but a physically organized circulation. The absence of the zonal extension in the two-dimensional framework intensifies the counter circulation.

Figure 3. Vertical-meridional cross section of temporal mean over the last 700 days of the 1,000 days time integration for the DY2-L128-DZ35 case. Contours are the mass-flux stream function at −2.0e9 (dashed) and 2.0e9 (solid). Smoothing is applied to two contour lines to improve readability.

Figure 4. Same figures as Figure 3, but for (a) vertical velocity, (b) mixing ratio of water vapor, (c) radiative heating rate, and (d) mixing ratio of total hydrometeors.

The simulated Hadley circulation in the two-dimensional framework has the same structure as the observed Hadley circulation seen in the reanalysis data in the sense of the zonal mean (Nguyen et al., 2013; Waliser et al., 1999). Air temperature (Figure 3d) is relatively homogeneous in the meridional direction from the Equator to 20°and has a gradient from 20° to 50°. This temperature profile is the result of a balance between temperature and the prescribed SST through surface fluxes, convection, and the Hadley circulation. A similar shape to the temperature profile exists in the mixing ratio of water vapor (Figure 4b). High temperature and high moisture in the Tropics drive deep convection where the concentrated upward motion and mixing ratio of hydrometeors are found, especially around 8° (Figures 4a and 4d). The updrafts associated with the deep convection and the poleward outflows in the upper troposphere transport a relatively higher equatorial angular momentum to the mid-latitudes (Figure 3b), which accelerates westerly winds (Hide, 1969; T. Schneider et al., 2010). The acceleration of the westerly wind reaches a maximum where the mass-flux stream function approaches zero around 40° (Figure 3c).

The pole-side edge of the simulated Hadley circulation is found around 45°, at the boundary of positive Ψ and negative Ψ (Figure 3a). This location is more poleward than in the real atmosphere; the pole-side edge of the real Hadley circulation is located from 30° to 35° according to various reanalysis data sets (Nguyen et al., 2013). In this simulation, a temporal mean zonal wind speed of $$ > $$60 ms⁻¹ is distributed from 30° to 45° (Figure 3c) and the maximum reaches 76 ms⁻¹. In the real atmosphere, a monthly mean wind speed for sub-tropical jet in September, which includes the autumn equinox, ranges from 40 ms⁻¹ to 50 ms⁻¹ over the region from 20°S to 30°S based on ERA-Interim data (Chenoli et al., 2016). Thus, in this simulation, the location of the zonal jet shifts poleward and has higher intensity than in the reanalysis data, which can be interpreted as a result of the absence of eddy diffusion due to no zonal waves in the vertical-meridional two-dimensional framework; the jet speed tends to be faster when the baroclinic waves are inactive (Nakamura, 1992). Although the locations and the intensity are different from those observed in the real atmosphere, we conclude that the two-dimensional framework produces a structurally equivalent Hadley circulation compared with the zonal mean of the observed real Hadley circulation.

In response to the Hadley circulation, various clouds are generated in the domain in addition to deep cumulus clouds. In the lower troposphere between 12° and 30° subtropical shallow clouds form in the subsidence-dominated environment, which suppresses deeper convective clouds (Figures 4a and 4d). These simulated shallow clouds are deeper than typical stratocumulus clouds and are mostly trade-wind cumulus. Typical stratocumulus is not well simulated, even in three-dimensional aqua-planet models due to the zonally symmetric boundary condition (Medeiros & Stevens, 2009). In the upper troposphere around 300 hPa height, there are weak signals of altocumulus or cirrus clouds in the hydrometeor field (Figure 4d). The vertical profile of radiative heating (Figure 4c) is different among these various clouds; a deep heating area exists in the area of the deep cumulus clouds and a shallow cooling region is found in the trade wind cumulus area above the shallow heating at the surface. Broad cooling is also found in the upper troposphere, which is associated with the anvil of the deep cumulus clouds or cirrus clouds. Thus, in the two-dimensional framework, the Hadley circulation generates various cloud types, associated with various physical processes.

Resolution Dependencies for Horizontal and Vertical Resolutions Sensitivities of the Simulated Hadley Circulation to Resolution

The characteristics of the simulated Hadley circulation explained in Section 3 are much affected by both horizontal and vertical resolution. Figure 5 shows the distribution of the mass-flux stream function for the different resolutions. The simulation case with the highest resolution among the nine (DY4-L128-DZ35; Figure 5a) is similar to the 2 km simulation of DY2-L128-DZ35 (Figure 3a) with respect to a pair of clockwise and counter-clockwise circulations over the Tropics. A simulation with high horizontal resolution tends to have an intense Hadley circulation for resolutions from 4 to 16 km; for example, the Hadley circulation of DY4-L128-DZ35 (Figure 5a) is stronger than that of DY8-L128-DZ35 (Figure 5d). The counter circulation at the Equator also becomes more intense at higher horizontal resolution at a given vertical resolution. Note that an opposite sensitivity is found when the horizontal resolution changes from 4 to 2 km; the circulation in DY2-L128-DZ35 (Figure 3a) is weaker than in DY4-L128-DZ35. A similar trend of strong counter circulation can be found for vertical resolution changes; in the DY4 cases, DY4-L128-DZ35 has a stronger counter circulation than DY4-L128-DZ50 and DY4-L128-DZ100. In the DY16 cases, the counter circulation is narrow and unclear, which is likely to be related to the broadening of the Hadley circulation to the equatorial area. These features can also be found in the hydrometeor field in Figure 6; at a given vertical resolution of L128-DZ35, the hydrometeor mixing ratio in DY16 tends to be less than in DY4, and the location of deep convection in DY16 shifts toward the Equator. The dependence of the Subtropical low level clouds on spatial resolution is not clear from Figure 6.

Figure 5. Vertical-meridional cross section of the mass-flux stream function for the L128 cases. Temporal means for the last 700 days of the 1,000 days time integration are plotted. The “⇔” symbol in the upper right corner means that the case is identified as a split convection case.

Figure 6. Vertical-meridional cross section of the mixing ratio of total hydrometeors for the L128 cases. Temporal means for the last 700 days of the 1,000 days time integration are plotted. Contours are the mass-flux stream function at −2.0e9 (dashed) and 2.0e9 (solid). Smoothing is applied to contours to improve readability. A “⇔” symbol at upper right corner means that the case is identified as a split-circulation case.

An interesting change in the simulated Hadley circulations is that the Hadley circulation is split into two parts around 15° at low resolution, especially obvious in DY8-L128-DZ100 and DY16-L128-DZ100 (Figures 5f and 5i, respectively). The split also appears in DY8-L128-DZ50 (Figure 5e). Hereafter, we refer to cases with the split as “split-circulation cases.” A “⇔” symbol in Figures 5 and 6 indicates the split-circulation cases, identified by the split index defined in Appendix C, which is based on the sign of the vertically integrated mass stream function. Similar characteristics are also observed in the hydrometeor field as shown in Figure 6. In low vertical resolution cases, the area corresponding to deep convection is split into two parts; for instance, deep convection is found around 7° and 18° in the DZ100 cases (Figures 6c, 6f, and 6i). DY4-L128-DZ100 (Figure 6c) and DY16-L128-DZ50 (Figure 6h) seem to be split cases in the hydrometeor field, while these are not identified as a split-circulation case by the split index. The split in the area of deep convection becomes clear at low horizontal resolution (e.g., Figures 6a, 6d, and 6g). It should be noted that the split is also observed in the case with high FA-resolution and low BL-resolution (DY4-L256-DZ100; not shown), which indicates a causal relationship between the split and the BL-resolution rather than with the FA-resolution. This will be discussed more deeply in Section 5.1.

In order to summarize the resolution dependencies for the simulated Hadley circulation, the relationships between the resolution, the pole-side edge of the simulated Hadley circulation, and the intensity of the westerly jet are diagnosed and plotted in Figure 7. Following Stachnik and Schumacher (2011) and Nguyen et al. (2013), the pole-side edge of the Hadley circulation is defined as the latitude with the first zero (changing from positive to negative) of the vertical summation of the mass-flux stream function from the Equator. The vertical summation is calculated from 800 hPa to 200 hPa to avoid contamination caused by conditions in the lower and upper troposphere. Then Gaussian smoothing is applied to the vertically accumulated mass-flux stream function to remove the small scale fluctuations in order to have a single zero around the pole-side edge of the Hadley circulation. While Nguyen et al. (2013) used a height range from 600 hPa to 400 hPa, we use a wider range in order to take into account the variability in the stream functions (Figure 5). The location of the pole-side edge of the simulated Hadley circulation is influenced by the height range of the vertical summation as in Table 3; comparing to the height range from 800 hPa to 200 hPa, the pole-side edge tends to reside at a higher latitude when vertical summation is computed with an upper height range from 500 hPa to 200 hPa and tends to be at a lower latitude when vertical summation is computed with a lower height range from 800 hPa to 500 hPa. In addition, the variation of the location of the pole-side edge tends to be small at high horizontal resolution. A westerly jet speed greater than 40 ms⁻¹ is used here as a proxy (i.e., index) for the intensity of the Hadley circulation in the two-dimensional simulation. As already noted, 40 ms⁻¹ is a typical value for a sub-tropical jet in the real atmosphere (Chenoli et al., 2016). Although the westerly jet in our ideal simulations differs from that in the real atmosphere (c.f., Figure 3c), we choose 40 ms⁻¹ as a threshold value to identify the area of the westerly jet. The grids with the westerly jet are identified over the height range from surface to 50 hPa height. In the figure, the values for the 25, 50, and 75 percentiles of the westerly jet speed are shown for each case.

Figure 7. The width of the simulated Hadley cell and westerly jet speed. For the westerly jet speed, the bottom whisker, the center marker, and the top whisker show 25, 50, and 75 percentiles, respectively. Markers and bars placed outside the panel show mean (markers) and distribution (bars) of latitude for each group in combinations of horizontal and vertical grid spacings.

Table 3 Location Sensitivity of the Pole-Side Edge of the Simulated Hadley Circulation to the Height Range Using a Vertical Summation of the Mass-Flux Stream Function

Cases	Height range [hPa]
	200–800	200–500	500–800
DY2-L128-DZ35	46.1°	46.7°	45.9°
DY4-L128-DZ35	45.9°	46.3°	45.8°
DY4-L128-DZ50	44.7°	45.9°	44.2°
DY4-L128-DZ100	45.5°	46.0°	44.9°
DY4-L256-DZ30	44.6°	45.3°	44.1°
DY4-L256-DZ50	46.2°	46.9°	46.0°
DY4-L256-DZ100	44.9°	45.8°	44.5°
DY8-L128-DZ35	46.4°	47.5°	43.5°
DY8-L128-DZ50	45.9°	48.4°	45.2°
DY8-L128-DZ100	46.0°	47.1°	45.5°
DY16-L128-DZ35	46.0°	48.7°	45.3°
DY16-L128-DZ50	44.4°	47.2°	43.5°
DY16-L128-DZ100	45.5°	47.6°	44.4°

Note. The values are latitude for the pole-side edges, and the height range of 200–800 hPa corresponds to the plots in Figure 7.

The high vertical BL-resolution seems to make the jet speed slower; it has the slowest median of the jet speed in Figure 7 except for the DY4-L256 cases. For example, the westerly jet speed in DY4-L128-DZ35 is slower than that in DY4-L128-DZ100. The increase in the FA-resolution also tends to reduce the westerly jet speed; the DY4-L256 has a slower speed than DY4-L128 except for the DZ100 case. For the pole-side edge of the simulated Hadley circulation, two weak trends are found from Figure 7. One of the trends is that the pole-side edge for the L128 cases converges at around 46° as the vertical BL-resolution increases (blue colored markers). This is a positive feature and motivation to use the two-dimensional framework for physics scheme development, since one can expect to have a reference solution with a high resolution simulation. The sensitivity to the vertical BL-resolution does not show clear convergence, while the distribution width of the DY4-L128 cases is narrower than those of the DY16-L128 cases (gray horizontal bars at the bottom of Figure 7).

For horizontal resolution dependency, the locations of the pole-side edge of the simulated Hadley circulation are affected by horizontal resolution, but no systematic trend is found in the two-dimensional simulations. Hohenegger et al. (2020) pointed out that the inter-tropical convergence zone in three-dimensional real atmosphere simulations tends to locate at a lower latitude in higher horizontal resolution simulations over the Atlantic ocean, while the trend is unclear over the Pacific ocean. In contrast, Wan et al. (2008) found that high horizontal resolution results in a poleward shift of the zonal jet in idealized three-dimensional model simulations. Overall, horizontal resolution dependency of the locations of the pole-side edge of the simulated Hadley circulation does exist; the trend is, however, unclear. This implies that aspects other than horizontal resolution may affect the location of the simulated Hadley circulation.

Sensitivities for Cloud Characteristics

We first define the Tropical region and the Subtropical region to investigate the cloud characteristics simulated in the two-dimensional framework. We do not investigate mid-latitude clouds due to the absence of zonal waves, which are an important driver of clouds in the region. The Tropical region is defined where a quasi-barotropic stratification is maintained. The Tropical region starts at 0° and ends at a grid that meets the condition of |∂T/∂φ| > ɛ. From 0°, |∂T/∂φ| is computed at 200 hPa height and the condition is checked with ɛ = 0.01. The Subtropical region is defined as the region between the Tropical region and the Mid-latitude region and the latter is defined as the latitudes where the zonal wind speed at 200 hPa height is faster than 40 ms⁻¹. Before defining these conditions, we took averages for temperature and zonal wind speed for the analysis period of 700-day, and for both hemispheres. The boundaries of the Tropics-Subtropics and the Subtropics-Mid-latitude identified by these definitions are summarized in Figure 8. The Tropics-Subtropics boundary ranges from 10.5° to 17.5°, and the Subtropics-Mid-latitude boundary ranges from 26.0° to 29.0°.

Figure 8. Locations of region boundaries for the simulated Hadley circulations. Closed markers show the boundary between the Tropics and the Subtropics, and open markers show the boundary between the Subtropics and the Mid-latitude. The vertical line corresponds to the Subtropical region.

We identify two trends in resolution dependencies from Figure 8. The first is that the Tropics-Subtropics and the Subtropics-Mid-latitude boundaries tend to be at lower latitude as horizontal resolution increases for a given vertical resolution, expect for the L128-DZ100 cases. The second trend is that the highest vertical BL-resolution tends to have the widest Subtropics region at a given horizontal resolution, except for the DY4-L256 cases.

To investigate cloud features in each region, statistics are computed for cloudy columns. The liquid cloudy column is defined with max(q_c(φ)) ≥ 0.01 g kg⁻¹, where q_c(φ) is the vertical profile of the cloud liquid water mixing ratio at φ. Clouds containing only ice particles (e.g., cirrus cloud) are excluded here because our focus is on low level clouds and/or convective clouds. The conditional statistics for liquid cloudy columns are computed for all grid columns over the 700-day analysis period in each region. To maintain consistency with other parameters, the hemispheric mean is also calculated.

For the liquid cloudy column fraction, opposite trends are found at a given vertical resolution and a given horizontal resolution. For horizontal resolution dependency at a given vertical resolution (Figures 9a and 9b), cloud occurrence tends to decrease at high horizontal resolution in both the Tropical and the Subtropical regions. DY8-L128-DZ35 in the Tropics region and DY2-L128-DZ35 and DY16-L128-DZ50 in the Subtropics region do not follow this trend. A similar trend of horizontal resolution dependency for simulated deep convection was reported in a three-dimensional global model (Kajikawa et al., 2016; Miyamoto et al., 2013); they found that deep convection tends to distribute sparsely and the area of clear sky tends to increase at high horizontal resolution. Prein et al. (2015) also reported that cloud cover in regional climate models without convective parameterization is smaller than that in conventional GCMs. Dueben et al. (2020) did not find a remarkable difference in a global model between a simulation with convective parameterization at 9 km horizontal resolution and a simulation without convective parameterization at 1.45 km of horizontal resolutions. Determining the effects of resolution is difficult based on a set of test cases with different configurations of the model.

Figure 9. Resolution dependencies for the liquid cloudy column fraction: (a), (b) are for horizontal resolution changes and (c), (d) are for vertical resolution changes. (a), (c) for the Tropics and (b), (d) for the Subtropics. Temporal mean for 700 days in the analysis period. A value of 1 means cloudy for all columns in the region and for the entire 700 days.

For vertical resolution dependency at a given horizontal resolution (Figures 9c and 9d), cloud occurrence tends to increase at high vertical BL-resolution both in the Tropics and the Subtropical regions; the liquid cloudy column fraction tends to be high at high vertical BL-resolution both in the Tropics region and in the Subtropics region, except for DY16-L128. In a conventional GCM, a similar trend of vertical resolution dependency for shallow convection has been reported both in the Tropics and the Subtropics (Bogenschutz et al., 2021; Lee et al., 2021). Xu and Cheng (2013) also found that the global mean low-cloud amount simulated by a GCM using super parameterization (i.e., a multiscale modeling framework; DeMott et al., 2010) increases at high vertical resolution and approaches the satellite-observed low-cloud amount. The counteracting effect on cloud occurrence for the horizontal and the vertical resolutions found in our two-dimensional simulations is likely to be found in future long-duration GHM simulations, since consistent characteristics are found in other studies that use three-dimensional models.

Figure 10 shows probability distributions at various resolutions for convective available potential energy (CAPE), surface precipitation, cloud water path including both liquid and ice, and cloud top height for the Tropics region. Figure 11 shows probability distributions in the Subtropics region in the same manner as Figure 10, but cloud water path in Figure 11c consists of only liquid cloud water in order to focus on the low level clouds. The cloud top height is defined at the height of the highest level at which a condensate path is greater than 0.01 kg m⁻²; the condensate path includes both liquid and ice for the Tropics in Figure 10d and includes only liquid for the Subtropics in Figure 11d. Since this study is the first to apply the two-dimensional framework to cloud system resolving modeling, we report findings in detail, so that they may serve as a reference for future studies. Since the discussion of these figures runs the risk of being long and tedious, we focus on the following themes: (a) horizontal resolution dependencies in the Tropical region, (b) vertical resolution dependencies in the Tropical region, (c) horizontal resolution dependencies in the Subtropical region, and (d) vertical resolution dependencies in the Subtropical region.

Figure 10. Boxplots in the Tropics region for (a) convective available potential energy, (b) surface precipitation, (c) cloud water path including both liquid and ice, and (d) cloud top height. Samples are for the 700-day analysis period in each region. (a) Includes all grids in the sample, but (b–d) include only the cloudy column. The upper label denotes horizontal grid spacing and number of vertical levels, and the bottom label along the abscissa axis denotes the vertical grid arrangement. A combination of the upper and bottom labels specifies one test case. The box indicates the 25 percentile (Q1) to 75 percentile (Q3) range for each sample. The black bar in each box is the median and the open black circle is the mean. The lower (upper) whisker shows the lowest datum above Q1 - 1.5*(Q3–Q1), and the upper whisker shows the highest datum below Q3 + 1.5*(Q3-Q1).

Figure 11. Same as Figure 10, but for the Subtropics. To focus on the low level clouds, the cloud liquid water path is shown for (c) and the cloud top height diagnosed with liquid cloud water is shown for (d).

Horizontal Resolution Dependencies in the Tropics

As seen in Figure 10a, the thermodynamical environmental conditions for convection tend to be less suitable at high horizontal resolution, which is understood by lower mean (open circle) and median (bar) CAPE at higher horizontal resolution at a given vertical resolution. The width of the distribution (between whiskers) also tends to be narrow at high horizontal resolution. This trend in CAPE is clear for the L128-DZ35 cases in particular, but it is weak for L128-DZ50 and L128-DZ100 cases; for the L128-DZ50 cases, mean and median in DY8-L128-DZ50 are lower than in DY4-L128-DZ50, and for the L128-DZ100 cases, the difference among the cases is small. Regarding surface precipitation, the width of the distribution becomes narrow and the value of the 75-percentile (top of the boxes) becomes small as the horizontal resolution increases (Figure 10b), so the fraction of weak precipitation increases at high horizontal resolution. The mean values for the cloud and ice water path (Figure 10c) and the cloud top height (Figure 10d) tend to be high at high horizontal resolution at a given vertical resolution, expect for L128-DZ50. This is consistent with what has been found in a three-dimensional GHM; Hohenegger et al. (2020) pointed out the occurrence of deeper convective cloud increases in high horizontal resolution simulations.

The trend of more weak precipitation at high horizontal resolution can be interpreted as an overestimate of precipitation amount at low horizontal resolution rather than an increase in weak precipitation occurrence at high horizontal resolution. The possible reasons for the trend of weak precipitation at high horizontal resolution are one or a combination of follows: lower frequency of convection, smaller area of convection, or weaker precipitation rate in convective regions. Identifying the exact mechanism is out of the scope of this study and postponed for future work. We note that the bias of larger precipitation at low horizontal resolution over the Tropics has been reported in three-dimensional simulations, especially over the open ocean (Hohenegger et al., 2020; Miyakawa & Miura, 2019; Wedi et al., 2020). Miyakawa and Miura (2019) suggested considering the rate of rainwater production; a high updraft velocity in convective cores at high resolution (3.5 km resolution) is one of the factors explaining the small amount of precipitation. Wedi et al. (2020) pointed out that a simulation at 9 km resolution without convective parameterization causes excessively strong deep convection in terms of convective gravity waves and the simulation at 1.4 km resolution provides a better representation of deep convection and storms. Hohenegger et al. (2020) also found a decreasing trend of 40-day mean precipitation over the Tropical ocean as horizontal resolution increases from 80 to 2.5 km. These trends reported by previous studies are consistent with the trend found in our two-dimensional simulations and provide confidence in the usefulness of the two-dimensional framework. However, Dueben et al. (2020) pointed out that a simulation at 9 km resolution with convective parameterization generated a 10%–15% lower precipitation amount than a simulation at 1.45 km resolution without convective parameterization. In a comparison between conventional GCMs and regional high resolution models without convective parameterization, a large difference was not found for daily mean precipitation, whereas intense precipitation was found in regional high resolution models at the sub-daily time scale (Prein et al., 2015). These results clearly tell us that the dependence of precipitation amount on resolution is affected by the convective parameterization.

Vertical Resolution Dependencies in the Tropics

For the vertical resolution dependencies in the Tropics, thermodynamic environmental conditions conducive to convection tend to occur at high vertical BL-resolution; the mean and median CAPE (Figure 10a) tend to increase as the vertical BL-resolution increases at a given horizontal resolution, except for the DY4-L128 cases. The mean and median of CAPE in the DY4-L256 cases is also higher than that in the DY4-L128 cases at a given vertical resolution (note, DY4-L256-DZ30 is compared to DY4-L128-DZ35). This implies that the high FA-resolution generates high CAPE. The trend observed in the liquid cloudy column fraction for vertical BL-resolutions (Figure 9c) is consistent with the increasing trend of CAPE if one assumes the high vertical BL-resolution increases the frequency of deep convection due to high CAPE. The frequency of deep clouds increases at high vertical BL-resolution; the 75-percentile value and the extreme value for the cloud and ice water path (Figure 10c) and the cloud top height (Figure 10d) become high at high vertical BL-resolution at a given horizontal resolution, except for DY4-L128-DZ50, DY4-L256-DZ100, and DY8-L128-DZ100. The 25-percentile value of the cloud top height also becomes high at high vertical Bl-resolution, except for DY8-L128-DZ50. However, the precipitation amount at convective region decreases at high vertical BL-resolution at a given horizontal resolution; the 75-percentile value for surface precipitation (Figure 10b) decreases as vertical BL-resolution increases. The high FA-resolution produces a systematic vertical resolution dependency. The mean and 75-percentile values of CAPE in DY4-256 increases, the mean of precipitation in DY4-L256 decreases, and the 75-percentile value of cloud top height increases systematically as vertical resolution increases.

Horizontal Resolution Dependencies in the Subtropics

Similarly to the horizontal resolution dependency in the Tropics, the thermodynamic environmental conditions for convection tend to be stable and the frequency of moist convection tend to decrease at high horizontal resolution: the mean (open circle) and median (bar) CAPE (Figure 11a) tend to be smaller at high horizontal resolution at a given vertical resolution. Since the trends in CAPE are hard to identify in the L128-DZ100 cases, they may become less susceptible at low vertical BL-resolution. The frequency of weak precipitation increases at high horizontal resolution; the 75-percentile of the surface precipitation (top of box in Figure 11b) clearly decreases as the horizontal resolution increases, while the mean surface precipitation does not seem to change (in L128-DZ35) or slightly increases (in L128-DZ50 and L128-DZ100) as horizontal resolution increases. The mean and 75-percentile value of the cloud water path (Figure 11c) also tend to decrease as the horizontal resolution improves, except for the L128-DZ100 cases. A similar decreasing trend is found for the 75-percentile value of the cloud top height for warm clouds (Figure 11d), but it is less obvious than that for the cloud water path. The tendency for less precipitation and lower cloud water path are reasonable under higher atmospheric stability, as indicated by less CAPE at high horizontal resolution, and are also consistent with the decreasing trend in low liquid cloud fraction (Figure 9b) at high horizontal resolution.

Vertical Resolution Dependencies in the Subtropics

In the Subtropics, the mean, median, and 75-percentile values of CAPE (Figure 11a) increase at high vertical BL-resolution at a given horizontal resolution, except for DY8-L128-DZ50, and this implies that the thermodynamic environmental conditions for convection tend to become more favorable at high vertical BL-resolution. However, mean and 75-percentile values of surface precipitation (Figure 11b) tend to be small at high vertical BL-resolution at a given horizontal resolution. The mean values of L128-DZ50 in DY8 and DY16 are exceptions; they have almost the same values as in L128-DZ100. In addition, mean cloud top height (Figure 11d) tends to be low at high vertical BL-resolution, except for the DY8-L128 cases. The vertical resolution dependency seems to be more marked in DY4-L128 than in DY8-L128 and DY16-L128; the decreasing trends with increasing vertical BL-resolution are found in the mean values for surface precipitation, cloud water path, and cloud top height for warm clouds in the DY4-L128 cases, while systematic trends are absent in the DY8-L128 and DY16-L128 cases.

Similarly to the Tropics, high FA-resolution makes CAPE even higher; mean, median, and 75-percentile values of CAPE in DY4-L256 are larger than those in DY4-L128 for comparable BL-vertical resolutions (comparison between L256-DZ30 and L128-DZ35, L256-DZ50 and L128-DZ50, and L256-DZ100 and L128-DZ100). While heavier total precipitation generally occurs under higher CAPE conditions due to more frequent convection, the mean surface precipitation in DY4-L256 is smaller than that in DY4-L128.

Summary for Boxplots

Considering the horizontal and vertical BL-resolution dependencies found in Figures 10 and 11, in response to spatial resolution improvements, in the Tropics, the frequency of deep convection increases (Figures 10c and 10d) and the intensity of convection, as represented by the surface precipitation, becomes weaker (Figure 10b). This trend of weaker precipitation at higher horizontal resolution is likely to exist in future long-duration GHM simulations, since many similar trends have been reported in three-dimensional models. In the Subtropics, high spatial resolution shifts the low clouds towards weaker and shallower convection (Figures 11b–11d). In our simulations, these Tropical and Subtropical responses to spatial resolution manifest as a narrower Tropical area at high resolution once a balanced state is achieved. To the best of our knowledge, a three-dimensional GHM simulation to support the resolution dependencies of low clouds in the Subtropics and the narrower tropical area at high resolution does not exist in literature.

For both regions, the convective instability decreases with increasing horizontal resolution (Figures 10a and 11a) while it increases at high vertical resolution, so that these effects counteract each other. In addition, we found that precipitation amount decreases as spatial resolution (horizontal and vertical) improves. These results are counterintuitive, as one might have expected high precipitation amounts with high convective frequency and instability. Unfortunately, we have not found literature sources to validate these results at this time, in particular for vertical resolution dependency. Further research is required to address the mechanisms behind the results and to determine a generality of this trend.

Discussion Why Does the Circulation Split at Low Resolution?

In Section 4, we showed the split in the simulated Hadley circulation for the cases with low resolution, particularly at low vertical resolution (Figure 5). Among all 13 simulations, the split in the simulated Hadley circulation is found in three cases (DY8-L128-DZ50, DY8-L128-DZ100, and DY16-L128-DZ100) and a weak split in the circulation is also seen in DY4-L128-DZ100, DY16-L128-DZ35, and DY16-L128-DZ50. It should be noted that a split in the simulated Hadley circulation is not found in the DY4-L256 cases.

Figures 12a and 12b show vertical-meridional cross-sections of temperature and water vapor mixing ratio comparing the split-circulation case of DY8-L128-DZ100 to the non-split case of DY8-L128-DZ35. In these figures, the split for DY8-L128-DZ100 is identified by contour lines of the stream function; an area with counter-clockwise circulation (dashed contour line) exists between 10° to 20°and is located between two clockwise circulations. The split case has two areas of deep convection at the upward branch of the two clockwise circulations at 6° and 18° (c.f., Figure 6f).

Figure 12. Vertical-meridional cross-sections of the temporal mean fields of (a) temperature and (b) water vapor mixing ratio; left is DY8-L128-DZ35 (non-split case) and right is DY8-L128-DZ100 (split case). Contours are the mass-flux stream function at −2.0e9 (dashed) and 2.0e9 (solid). (c) The temporal mean of CAPE for several cases. The temporal mean is taken over the last 700 days of the 1000-day time integration.

In the split-circulation case, deep convection is not likely to occur in the area where the downward branch of the clockwise circulation and that of the counter clockwise circulation meet (e.g., approximately 12° for DY8-L128-DZ100), that is., where the mean vertical motion is expected to be downward and a convectively stable condition is anticipated. For DY8-L128-DZ100, the column mean vertical speed in the level from 850 hPa to 200 hPa is indeed downward in this area. Around 12° latitude, the temperature (Figure 12a) below 900 hPa height in the split-circulation case (DY8-L128-DZ100) is more homogeneous and warmer than that in the non-split case (DY8-L128-DZ35) and the water vapor mixing ratio (Figure 12b) in the split-circulation case is less than that in the non-split case at 12°. These temperature and water vapor characteristics are common among the split-circulation cases at other horizontal and vertical resolutions. This warmer and drier conditions in the lower troposphere are likely generated by subsidence associated with the large-scale circulation, which suppresses moist convection by establishing relatively stable conditions. These are the conditions associated with the split-circulation case. A cause of the subsidence could be the appearance of two deep convective regions at 6° and 18° (c.f., Figures 5f and 6f) that organize the subsidence between them. The thermodynamic conditions that organize these two deep convective regions at 6° and 18° are found in the meridional profile of CAPE (Figure 12c) for the split-circulation cases (DY8-L128-DZ100 and DY16-L128-DZ100); one is the maximum of CAPE around 7°corresponding to the deep convection at 6°and the second is the local peak in CAPE around 17° corresponding to the deep convection at 18°. On the other hand, CAPE decreases from around 6°–20° monotonically for the non-split cases (DY4-L128-DZ35 and DY8-L128-DZ35). There is an exception in the weak split case; while both DY4-L128-DZ100 and DY16-L128-DZ35 are the weak split cases, DY4-L128-DZ100 does not have a local peak.

To clarify the cause of the split circulation from a different perspective, we calculate power spectra for u (zonal wind speed) at 200 hPa height, v (meridional wind speed) at 200 hPa and 850 hPa heights, and w (vertical velocity) at 850 hPa height. Figure 13 shows the meridional distributions for these power spectra for u and v at 200 hPa height for the DY4, DY8, and DY16 cases and Figure 14 shows the same figures but for v and w at 850 hPa height.

Figure 13. Meridional distributions of the power spectra at 200 hPa height for u (contour) and v (shaded). The power spectrum is calculated using data for 700 days in the analysis period at each meridional grid point, and averaged over both hemispheres.

Figure 14. Meridional distributions of the power spectra at 850 hPa height for w (contour) and v (shaded). The power spectrum is calculated using data for 700 days in the analysis period at each meridional grid point and averaged over both hemispheres.

A peak in the power spectrum of u (line contours in Figure 13) is found at about 35° at a frequency >$$1{0}^{-2}$$ day⁻¹ (i.e., 100-day period). Considering the latitudinal location and the height, the peak in u is associated with the westerly jet corresponding to the simulated Hadley circulation. This fluctuation in u is related to the global scale, which is indicated by a fluctuation with a 100-day period found in the temporal evolution of the global accumulated angular momentum in Figure 2. The power spectrum of u has clear resolution dependency both for horizontal and vertical resolutions; the power spectrum of u tends to expand to lower latitudes at low horizontal resolution and at low vertical BL-resolution. The trend in the vertical BL-resolution in the DY4 cases is not as significant as that in the other horizontal resolution cases.

The power spectrum of v at 200 hPa height (shaded contours in Figure 13) shows a remarkable difference between the split-circulation cases (indicated by a “ ⇔” symbol) and the non-split cases; the peak at high frequency around 2 × 10⁻¹ day⁻¹ (i.e., 5-day period) tends to be weak and the peaks in the low frequency around 10⁻² day⁻¹ tend to be strong in the split-circulation cases, while the high-frequency peak tends to be strong and the low-frequency peaks tend to be very weak in the non-split cases. The reason has not been identified, but the high-frequency peak is quite weak in DY16-L128-DZ35. The high-frequency peak probably corresponds to an organized system with deep convection, which is indicated by the location between 0° and 15° with a 3–5 day period. The low-frequency peak likely corresponds to the simulated Hadley circulation, based on the overlap with the peak in u.

Like v at 200 hPa height, both the high-frequency peaks (2 × 10⁻¹ day⁻¹) and the low-frequency peaks (around 10⁻² day⁻¹) are found for v at 850 hPa height (shaded contours in Figure 14). These low frequency peaks and their vertical BL-resolution dependence seem to be closely related to the split in the simulated Hadley circulation, as evidenced by the interesting low-frequency peaks at about 18° where the split circulation is located. The peaks at low frequency are found in all nine cases; they are, however, limited to the Equator (0°–5°) in the high vertical BL-resolution cases (Figures 14a, 14d, and 14g) and concentrated over two latitudinal locations in the equatorial area from 0° to 10° and in the off-equatorial area around 18° in the low vertical BL-resolution cases (Figures 14f and 14i). A horizontal resolution dependence is also found for the low-frequency peaks; the low frequency peaks tend to be weak at high horizontal resolution (e.g., Figure 14a). The horizontal resolution dependence implies that fluctuations over a short time period tend to become dominant in the high horizontal resolution cases.

Low-frequency peaks (10⁻² day⁻¹) around 18° are also found in the power spectrum of w at 850 hPa height (line contours in Figure 14) for the split-circulation cases (Figures 14e, 14f, and 14i). The wide range in temporal scales, from the single day period to the hundred-day period, is a common characteristic among all cases in w at 850 hPa height, however, the latitudinal locations of peaks can be categorized into two types; 1) the concentrated peaks over an area from 5° to 10° and 2) the split peaks located over the two latitudinal locations around 6° and around 18°. The former type is seen in the non-split cases (e.g., Figures 14a, 14b, and 14d) and the latter type is found in the split-circulation cases. The latter type can be also found in DY4-L128-DZ100 and DY16-L128-DZ35, which are considered weak split cases. Another trend is that the peaks in w at 850 hPa height tend to expand toward the Equator in the DY16 cases and this trend could be related to the weak counter circulation in the DY16-L128 cases shown in Figures 5g–5i.

We hypothesize that the split in the simulated Hadley circulation is fundamentally related to low-frequency (large-scale) processes rather than to high-frequency (convective) processes. One reason is that major resolution dependencies for the power spectra of v are found at a low frequency. Furthermore, overlaps of the low-frequency peaks (100-day period) between v and w at 850 hPa height are found at the latitudinal location of 18° where the counter clockwise circulation is found in the split-circulation cases (DY8-L128-DZ50, DY8-L128-DZ100, and DY16-L128-DZ100). This coincidence implies a particular interaction between v and vertical motion through low-frequency processes in the low vertical BL-resolution cases. In addition, the interaction between v and u through the Hadley circulation should be taken into account; for example, overlaps in the peaks between v and u at 200 hPa height are found at the low frequency in the split-circulation cases, particularly in DY16-L128-DZ100. Therefore, the non-split and the split simulated Hadley circulations is related to different balanced states via different scale interactions.

The reason why a combination of low horizontal resolution and low vertical resolution (e.g., DY16-L128-DZ100) has a clear split is explained by resolution dependencies as follows; as discussed in Section 4, low horizontal resolution as well as low vertical BL-resolution tend to generate intense deep convection, which leads to the formation of a subsidence region once two convective regions start to emerge. At the same time, the location of the deep convection areas become fixed, as shown by the spectral analysis, which results in a CAPE profile with two peaks. A linear response among the cases with a nearly fixed grid aspect ratio (e.g., DY4-L128-DZ35 (13a), DY8-L128-DZ50 (13e), and DY16-L128-DZ100 (13i)) can be observed, for example, the weakening trend of the low-frequency peaks and the intensifying trend of the high-frequency peak for v (Figure 13). This indicates the impact of the grid aspect ratio on the simulated Hadley circulation.

The split bias of the simulated Hadley circulation found in this study may not occur in three-dimensional simulations due to zonal variability in upsidence and subsidence, which act to avoid the locking of the location of intense subsidence. However, the difference in the scale interactions corresponding to the resolution dependencies found in the two-dimensional simulation can also be found in three-dimensional simulations. Following the hypothesized scenario, the well known bias of the double inter-tropical convergence zone pattern of precipitation in three-dimensional GCMs (Kodama et al., 2015; Lin, 2007; Tian & Dong, 2020) could be affected partially by variability in scale interactions between convection and circulation due to sensitivity to both horizontal and vertical resolution. In this study, we found an interesting resolution dependency such that clouds in the Subtropical region tend to be intense and relatively deep, and the Tropical region tends to expand to a relatively high latitude when the horizontal and vertical resolutions are low.

Impacts Associated With the Thickness of the Lowest Layer

We briefly mention the simulation cases with vertical grid arrangement of L128-DZ50z1 and L128-DZ100z1, which are aimed at assessing the impacts of the thickness of the lowest layer. We find that in the L128-DZ50z1 and L128-DZ100z1 cases there are no notable difference in trends when comparing to the vertical grid arrangement of L128-DZ50 and L128-DZ100. The thickness of the bottom layer has an impact on the location of the pole-side edge of the simulated Hadley circulation, but no systematic differences in L128-DZ50z1 and L128-DZ100z1 are found when comparing to L128-DZ50 and L128-DZ100. As in the vertical grid arrangement cases for L128-DZ50 and L128-DZ100, a split in the simulated Hadley circulation is found for DY4-L128-DZ100z1, DY8-L128-DZ50z1, DY8-L128-DZ100z1, DY16-L128-DZ50z1, and DY16-L128-DZ100z1. For cloud properties, the larger thickness of the lowest layer tends to generate a lower cloud fraction in both the Tropics and Subtropics and they tend to precipitate less in the Tropics compared to the vertical grid arrangement cases for L128-DZ50 and L128-DZ100. The smaller thickness of the lowest layer may prevent the moisture mixing up efficiently in the vertical direction or not allow the surface processes to generate strong enough moisture flux. The causes of these resolution dependencies are not clear at this time and they may differ between numerical models and/or in combination with other physics schemes and experimental settings.

Summary

We have proposed a framework for global model development geared toward extremely high resolution simulations. High costs of computation and large data volumes are inexorable challenges that burden the development and tuning of physics schemes in high resolution global models. Our solution to this problem is to remove the zonal extension of the global simulation, using a vertical-meridional two-dimensional simulation over a global pole-to-pole computational domain. The two-dimensional framework significantly reduces computation time and output data volume. It simulates a large-scale circulation corresponding to the Hadley circulation, which is a key feature of the Earth's climate. Importantly, by using fine grid spacing and by leveraging freed computational costs created by eliminating the zonal extension, we can explicitly (without a cumulus parameterization) simulate various cloud types embedded in the simulated Hadley circulation. This is an advantage of the two-dimensional framework since interaction between clouds and large-scale circulation is not attainable with single column and regional simulations. Such a two-dimensional framework could also be applied to ocean modeling.

We justify our approach by providing examples of some common behaviors between our two-dimensional simulation and three-dimensional simulations from previous studies; these similarities can be considered assessment points for physics scheme development. Furthermore, the differences between two- and three-dimensional simulations may help to elucidate processes in three-dimensional simulations.

We have demonstrated the applicability of this framework by conducting a study of the sensitivity of Hadley circulation and clouds to spatial resolution. A series of 1000-day simulations were performed with various horizontal and vertical grid spacing ranging from 16 to 2 km in horizontal grid spacing and from 35 to 100 m in the vertical grid spacing below 2 km height with 128 or 256 vertical levels. Some of the characteristics of the simulated Hadley circulation converge at high vertical resolution; the pole-side edge of the simulated Hadley circulation seems to converge at around 46° in latitude and the westerly jet speed tends to become slower at the highest vertical resolution. The trends for horizontal resolution dependency both in the pole-side edge and the westerly jet speed are not clear.

An interesting result found in this study is the split in the simulated Hadley circulation at low resolution. We found that the low resolution forcefully fixes the location of two deep convective areas and generates a counter circulation with subsidence in between. A power spectrum analysis of the wind speeds reveals a resolution dependence of the scale interactions; the Hadley circulation is governed by a high-frequency processes such as convection in the simulations with high resolution, while it is regulated by low-frequency processes in the simulations with low resolution. The low-frequency fluctuation in the wind fields is likely related to the split in the deep convection. With this new understanding, we further hypothesize that the well known bias of the double inter-tropical convergence zone pattern of precipitation in three-dimensional GCMs/GHMs is partially affected by the enhanced scale interactions at low frequency modes due to use of low horizontal and vertical resolutions. Note that the bias of the double inter-tropical convergence zone pattern has also been observed in a GHM (Kodama et al., 2015; Satoh et al., 2019).

We have found several resolution dependencies that are consistent with those in three-dimensional GHMs and we expect these trends to be found in future long-duration GHM simulations. The decreasing trend of precipitation amount as horizontal resolution increases is consistent with three-dimensional GHM results (Hohenegger et al., 2020; Miyakawa & Miura, 2019). The trend of an increase in deeper convection as horizontal resolution increases is also similar to that found in the three-dimensional GHM study by Hohenegger et al. (2020). Regarding the resolution dependence of the simulated cloud properties, both for the Tropical and Subtropical regions, the fraction of liquid cloudy columns tends to decrease at high horizontal resolution, while it tends to increase at high vertical resolution. This counteracting trend might be found in future long duration GHM simulations since the decreasing trend of cloud amount at high horizontal resolution and the increasing trend of cloud amount at high vertical resolution have been reported in previous three-dimensional studies. To comprehensively validate this counteracting trend, we need to test both horizontal and vertical resolution sensitivities in a GHM. These trends in cloudiness are related to the resolution dependence of environmental conditions; mean CAPE tends to decrease as horizontal resolution increases and increases as vertical resolution increases both in the Tropical and Subtropical regions. Considering horizontal and vertical resolution dependencies together, clouds and convection in the Subtropical region become weaker and shallower with increasing horizontal and vertical resolution. In addition, the mean cloud top height in the Subtropical region tends to be low at both high horizontal and high vertical resolutions. These trends in cloud properties provide important guidance for adequately representing shallow clouds in future GHM climate simulations, but these have not been validated by other models at this time.

The largest disadvantage of the two-dimensional simulation framework is the lack of zonal fluctuations like baroclinic waves and Rossby waves. This prevents one from studying the influence of eddy diffusion on the Hadley circulation, which has an important impact on the width of Hadley circulation in the real atmosphere. Clouds and convection are also affected by these synoptic scale waves and the resolution dependency could be more complicated in the three dimensional simulation. We need to wait for future studies that examine long duration GHM simulations and/or vertical-meridional two-dimensional simulations by other models to understand how horizontal waves and eddies affect the resolution dependencies found in this study. Even with these caveats, the framework can work well as a test-bed for performance assessments, and provide clues to understand the behaviors of physics schemes in global high resolution atmospheric simulations.

A AppendixImplementation of SHOC With Explicit Diffusion

SHOC is a non-local subgrid scale mixing scheme, which simulates diffusion assuming a single overturning event due to a dry or moist convective process. For each column, SHOC identifies cloud-free layers for dry convection and cloudy layers for moist convection, where each layer consists of one or more levels. The overturning process in each layer has a characteristic time scale (τ), a diagnostic variable used to determine a mixing length scale (L), where L ∝ τ^1/2, in a very simplified manner. Another time scale, a return-to-isotropy timescale (τ_v), is also diagnosed and used in the calculation of the eddy viscosity K = P_rK_h where P_r is the Prandtl number (P_r = 1 in SAM), K_h is the eddy thermal viscosity defined as K_h = C_kτ_ve, C_k = 0.1, and e is turbulence kinetic energy. An advantage of SHOC is that it only predicts e based on the 1.5-order TKE closure scheme (TKE1.5; Deardorff, 1980) with [Image Omitted. See PDF]where B, D, and S are terms associated with buoyant production, dissipation, and shear production, respectively, Δt is the time integration interval, and n is the timestep. B and D depend on L and are defined as [Image Omitted. See PDF]and [Image Omitted. See PDF]where θ_v is virtual potential temperature, w is vertical velocity, and $$\bar{{w}^{\prime }{\theta }^{\prime }}$$ is virtual potential temperature flux. We see that B ∝ τ_v when ∂θ_v/∂z < 0 and D ∝ τ^−1/2 since K_h ∝ τ_v and L ∝ τ^1/2. With this formulation, SHOC can generate larger e if one assigns reasonable time scales for a single convection event within a layer by making B larger and D smaller (and thus larger K compared with TKE1.5). We find that these time scales diagnosed by SHOC are typically on the order of 500 s. This is the solution that SHOC uses to address a known underestimation of mixing due to small e and K for TKE1.5 (Bogenschutz & Krueger, 2013; Cheng et al., 2010).

However a large K can be problematic when SHOC is coupled with an explicit time integration with a 10 s timestep since it potentially gives an unstable solution through over-mixing. This is the problem we encountered during implementation of SHOC into SAM. By performing thorough experimentation with various approaches, we arrived at our solution to the overestimated K for explicit diffusion.

Our solution to the problem is two-fold: (a) removal of τ_v from K; and (b) limiting e by re-scaling L. Instead of the definition of K used by SHOC, we use the definition originally proposed by TKE1.5: [Image Omitted. See PDF]where L* is the re-scaled length scale. Since the maximum time scale for mixing allowed in an explicit time integration is Δt, L diagnosed by SHOC is re-scaled as [Image Omitted. See PDF]when Δt > τ, we set L* = L. We also further limit L* by the condition l ≤ L* ≤ 2Δz, where l is the mixing length given by TKE1.5. This is the most relaxed setting and is practically stable for all cases in our two-dimensional Hadley simulations. With the re-scaled, limited L*, the buoyant production and dissipation terms are computed as [Image Omitted. See PDF]and [Image Omitted. See PDF]when Δt < τ, B* < B and D* > D are achieved since $${B}^{\ast }\sim B({\Delta}t/\tau )$$ and $${D}^{\ast }\sim D(\tau /{\Delta}t)$$ and e and K are commensurately smaller. At the lower limit of L* (i.e., L* = l), the modification makes SHOC behave like TKE1.5. On the other hand, for sufficiently large Δt and Δz, the modified SHOC becomes closer to the original SHOC and the differences between the original SHOC and the modified SHOC come solely from the difference in the definition of K.

B AppendixPerformance of the Modified SHOC

The performance of the modified SHOC is tested by simulating a trade wind cumulus case based on the Barbados Oceanographic and Meteorological Experiment (BOMEX; Siebesma et al., 2003) in a three-dimensional domain. Three types of turbulence schemes are compared; the modified SHOC, the original SHOC, and TKE1.5. Two horizontal grid spacings of 2 and 8 km are used with 64 horizontal grids (i.e., 128 × 128 km² and 512 × 512 km² domains). Some settings in this test are different from the configurations specified by the intercomparison project of the Global Energy and Water Cycle Experiment (GEWEX) Cloud Systems Study (GCSS); interactive radiation, interactive surface flux, and precipitation are calculated, a deeper domain (the L128-DZ100 vertical grid defined in this study with a domain top at 35 km) is used, and the prescribed large-scale forcing is not used. P3 and RRTMG are used to represent microphysics and radiation processes, respectively. Aerosol number concentration is fixed at 100 mg⁻¹. We do not use the diagnosed SGS cloud fraction and cloud mass from SHOC for this test. The duration of the simulation is 10 days with a 30 s time step except for the original SHOC, which uses a 4 s time step due to the numerical instability discussed in Appendix A. As a reference, a large-eddy simulation (LES) with 100 m horizontal grid spacing and a 6.4 × 6.4 km² domain is performed. The time step for the LES is 4 s. Additionally, another LES with the modified SHOC is carried out with 100 m horizontal grid spacing in order to examine the convergence of the modified SHOC to TKE1.5 for fine resolution and a small time step.

Figure B1 shows vertical profiles averaged over the computational domain for the last 7 days of the simulation period. The performance of the modified SHOC exhibits both improvements and degradations compared to the original SHOC and the LES and the performance places it between the original SHOC and TKE1.5 (except for the case of the modified SHOC with 100 m grid spacing). The turbulence kinetic energy below 3 km height in the original SHOC is closer to the LES than that in the modified SHOC. On the other hand, for the profile of turbulence kinetic energy above 3 km height, the performance of the modified-SHOC is better than or equivalent to the original SHOC. Similar trends are found in precipitation flux (Figure B1c); the performance of the modified SHOC above 4 km height is similar to the original SHOC (it is better with 8 km grid spacing), while the performance below 4 km height is degraded from the original SHOC. Furthermore, the numerical solution of the modified SHOC converges to TKE1.5 as grid spacing becomes finer, which is expected based on the design of the modified SHOC; the LES of the modified SHOC (mod-S 100 m in Figures B1a–B1d) is the closest to the LES.

Figure B1. Vertical profiles of temporal and domain mean values for the Barbados Oceanographic and Meteorological Experiment test case in a three-dimensional domain. The temporal mean was calculated over the last 7 days of the 10 days time integration. In the legend, turbulence kinetic energy (TKE), SHOC (Simplified High Order Closure), and modS mean TKE 1.5-order closure scheme, the original SHOC, and the modified-SHOC, respectively.

We found two concerns with the original SHOC through this BOMEX test. One is the differences in the buoyant production from the LES (Figure B1b). In the cloud layer between 1 and 3 km, the sign is opposite for two tested horizontal grid spacings and there are large positive and negative fluctuations above 3 km for 8 km horizontal grid spacing (Figure B1b). The turbulence kinetic energy is significantly larger than the LES above 3 km height for the original SHOC with 8 km horizontal grid spacing in (Figure B1a). There is no such bias in the modified SHOC with 8 km horizontal grid spacing.

The other concern is that the diagnosed SGS cloud fraction has a significant oscillation with time (not shown) both in the original SHOC and the modified SHOC when a pseudo single column model simulation is performed by copying horizontally one grid point to the other grid points both in x and y directions after each process calculation using the minimum grid number (4 × 4) that SAM allows. For this pseudo single column model test, the temporal mean values of the SGS cloud fraction in the original SHOC and the modified SHOC are at approximately 0.14. Even though SHOC can diagnose a SGS cloud fraction and SGS cloud mass, these SGS cloud properties should have a near constant value with time when an initially horizontally homogeneous environmental condition like BOMEX is prescribed. The oscillation is also found in the original code where the SGS cloud mass is coupled with a microphysics scheme and using Pr = 0.1 to compute eddy viscosity (Bogenschutz & Krueger, 2013). The version of Bogenschutz and Krueger (2013) blows up for the two-dimensional Hadley circulation simulations. Cheng et al. (2010) reported a similar oscillation in their third-order turbulence closure model, which they solved by parameterizing liquid water buoyancy. In this study, we do not use the SGS clouds diagnosed by SHOC to simplify implementation.

C AppendixDefinition of the Split Index

A split index (SI) is proposed to determine objectively the split in the simulated Hadley circulation instead of relying on human identification. A simulation case that meets the condition of SI > 0 is identified as a split-circulation case. SI consists of the summations of vertically accumulated mass-flux stream function as below, [Image Omitted. See PDF]where ∑_j and ∑_k are the meridional summation and the vertical summation, respectively; j is the meridional grid index and k is the vertical grid index. The height range for the vertical summation is from 800 hPa to 200 hPa; this height range is the same as that used in the analysis of the pole-side edge of the simulated Hadley circulation. The latitudinal range of the meridional summation for both the numerator and the denominator are taken over the area of the simulated Hadley circulation, but only negative values of the stream function are accumulated in the numerator.

D AppendixResolution Dependencies for the Simulations With the TKE1.5 Scheme

The two-dimensional simulations using the 1.5-order TKE closure scheme are presented for the 9 cases of the combinations of the horizontal grid spacings of DY4, DY8, and DY16 and the vertical grids of L128-DZ35, L128-DZ50, and L128-DZ100. The simulation settings are the same as the simulations in this study, except that we use the 1.5-order TKE closure scheme instead of the modified SHOC.

Figure D1, Figure D2 and Figure D3.

Figure D1. Same as Figure 7, but for simulation cases with the 1.5-order TKE closure scheme.

Figure D2. Same as Figure 8, but for simulation cases with the 1.5-order TKE closure scheme.

Figure D3. Same as Figure 9, but for simulation cases with the 1.5-order TKE closure scheme.

Acknowledgments

The authors thank to Dr. David A. Randall for his critical comments. Thanks also to the members of the Computational Climate Science Research Team, RIKEN Center for Computational Science, and Research Center for Urban Safety and Security, Kobe University, for their input. This work was funded by Scientific Discovery through Advanced Computing (SciDAC), award number DE-SC0018650, by the U.S. Department of Energy office of Biological and Environment Research. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. The authors acknowledge the NOAA Research and Development High Performance Computing Program for providing computing and storage resources that have contributed to the research results reported within this paper. ERA5 data set was generated and downloaded using Copernicus Climate Change Service information (2019).

Data Availability Statement

Data of this study are available at https://csl.noaa.gov/groups/csl9/datasets/data/cloud_phys/.