Reproducibility of Equatorial Kelvin Waves in a Superparameterized MIROC: 1. Implementation and Verification of Blockwise‐Coupled SP‐MIROC

The potential scope of superparameterization (SP) was extended to higher resolutions of the global climate model (GCM) component by devising a technique called blockwise coupling. In this method, a horizontal average of multiple GCM columns, instead of one, is coupled to a cloud‐resolving model (CRM) domain. This enables SP‐GCMs to reduce the computational cost drastically, enabling higher‐resolution GCMs to be superparameterized. A blockwise‐coupled SP‐GCM called SP‐MIROC was implemented by coupling the climate model MIROC6 to the CRM SCALE‐RM. The 4 × 4‐bundled SP‐MIROC successfully reproduced horizontal patterns and frequency distributions of precipitation and realistic amplitudes of equatorial Kelvin waves (EKWs), which were underestimated in the standard MIROC6. As discussed in Yamazaki and Miura (2024b, https://doi.org/10.1029/2023MS003837) of this study, the amplitude boost of EKWs was enabled by a top‐heavy heating in SP‐MIROC. Comparison of power spectra between the 4 × 4‐bundled SP‐MIROC and nonbundled SP‐MIROC indicated that the effective resolution of dynamic variables was not degraded by the blockwise technique. Rather, spectra in the 4 × 4‐bundled SP‐MIROC were more realistic than those in the nonbundled SP‐MIROC. Although the 4 × 4‐bundling limits convective coupling in the smallest GCM scale, it could offer the best match of resolutions between the GCM‐handled dynamics and SP‐derived physics because the effective resolution of the dynamics is lower than the nominal grid spacing.


Introduction
Atmospheric phenomena have wide range of spatial scales from microns to the global scale, and all of these phenomena are connected through various interscale processes.Therefore, accurate representation of global and synoptic features requires the appropriate handling of small-scale activities, such as cumulus convection and turbulence.
To incorporate the impact of cumulus convection on large-scale circulation, the traditional global climate models (GCMs) have deployed parameterization schemes.For example, Arakawa and Schubert (1974) developed a cumulus parameterization scheme by assuming a state of quasi-equilibrium between large-scale forcing and environmental feedbacks by cumulus plumes.Cumulus clouds were assumed to modify the environment through cloud top detrainment and inter-cloud downdrafts compensating for intracloud updrafts.They characterized subgrid cumulus clouds by normalized entrainment rates, which regulate cloud top heights and mass flux profiles.This assumption, combined with the aforementioned quasi-equilibrium closure, allowed them to diagnose forcings of subgrid convection from grid-scale environmental fields.However, the ability of such schemes to diagnose cumulus activities is not perfect.In the tropics, where subgrid clouds have a particularly large influence on large-scale circulations, this imperfection leads to difficulties in reproducing prominent large-scale disturbances (Bartana et al., 2022;J.-L. Lin et al., 2006) such as the Madden-Julian Oscillation (MJO) and equatorial waves.Furthermore, most GCMs have a bias pattern in the surface precipitation commonly called "double-ITCZ," in which the intertropical convergence zone (ITCZ) in the south-eastern Pacific is excessively strong (J.-L.Lin, 2007).
These issues motivated the development of global cloud-resolving models (GCRMs) or global convectionpermitting models, such as NICAM (Satoh et al., 2008;Tomita et al., 2005), ICON (Zängl et al., 2015), and MPAS (Skamarock et al., 2012); in these models, cumulus drafts and associated microphysical processes are explicitly treated in relatively high resolution.Many more models have been developed and joined intercomparison projects, such as DYAMOND (Stevens et al., 2019).Although supercomputers today enable ambitious modeling targeting climate simulations at 1-km resolution (Voosen, 2020), such models require large computational resources, making them unsuitable for long-term experiments.Superparameterization, originally termed Cloud Resolving Convection Parameterization (Grabowski, 2001;Grabowski & Smolarkiewicz, 1999), is a technique that incorporates the influences of cumulus convections explicitly with less computational costs than GCRMs.In a superparameterized (SP) model, a domain of a highresolution, cloud-resolving model (CRM) is embedded in every column of a GCM, and this embedded domain handles small-scale processes, including cumulus convections, while most of the atmospheric parameterizations are turned off in the GCM component.To reduce computational costs, the embedded CRM domains are usually two-dimensional and/or considerably smaller than the horizontal extent of the corresponding GCM grid.Even though two-dimensional CRMs behave differently to full 3D models, they feedback to the large scale similarly (Grabowski et al., 1998), except for high-frequency variabilities relatively irrelevant to the GCM dynamics.Through explicit representation of cumulonimbi, superparameterization is expected to enhance the reliability of climate simulations (Randall et al., 2003), although there remain problems, such as shallow clouds and boundary layer dynamics, in superparameterized models as well.It could also serve as a convenient framework to analyze multiscale interactions between the synoptic environment and cumulus convections because it can reproduce convection-dependent phenomena reasonably and handle mutual forcings between cumulus convection and large-scale environment explicitly.In contrast, GCRMs handle a wide range of scales seamlessly, making it difficult to draw distinct scale boundaries for analyses.Superparameterization belongs to the group of multiscale modeling techniques including the multiscale asymptotic method, in which differential operators and variables are expanded by a parameter representing the ratio between two scales of interest (Grooms & Julien, 2018).While the present superparameterization is not based on the multiscale asymptotic method and is sometimes described as an ad hoc coupling framework (Grooms & Julien, 2018), its integration to the multiscale asymptotic method can deepen the mathematical foundation to couple large and small scales more consistently.Superparameterization has been implemented to multiple GCMs including SP-CAM (Khairoutdinov & Randall, 2001), SP-CCSM (Stan et al., 2010), SP-EMAC (Jöckel et al., 2010), SP-WRF (Tulich, 2015), and E3SM-MMF (Hannah et al., 2020).It can reproduce tropical phenomena such as MJOs (Benedict & Randall, 2011;Hannah et al., 2015Hannah et al., , 2020;;Khairoutdinov et al., 2005;Stan et al., 2010) and equatorial Kelvin waves (EKWs) (Benedict & Randall, 2011;Hannah et al., 2020) well.Furthermore, it outperforms the conventional GCMs in reproducing mesoscale convective systems (MCSs, Kooperman et al., 2013), especially in high-resolution GCM configurations (G.Lin et al., 2022).Superparameterization can sometimes expose inhomogeneities in the GCM component as unrealistic convective patterns.Hannah et al. (2020) identified a spurious pattern called "grid imprinting" in the mean precipitation fields simulated by E3SM-MMF.This pattern consisted of lines of enhanced precipitation arranged in grids, reflecting the inhomogeneity of the spectral element grid used for dynamics.According to water budget analyses, the artificial circulation accompanying the grid imprinting had a non-negligible impact on the global mean.This grid imprinting bias, with smaller amplitude than that in E3SM-MMF, was also detected in the standard E3SM.Hannah et al. (2021) tackled against the grid imprinting bias in the standard E3SM by adopting a physics grid different from the dynamics grid.The E3SM dynamics employ the Galerkin spectral finite element method (Dennis et al., 2012).In this method, spectral elements containing 4 × 4 node points fill the globe.The node points, which are defined as quadrature points of polynomials consisting the spectral elements, are irregularly spaced and used as locations where the physics processes are evaluated.The irregularity of the physics grids leads to horizontally discontinuous physics tendencies, resulting in the grid imprinting bias.To solve this problem, Hannah et al. (2021) performed physical parameterization on a new finite-volume physics grid featuring quasiequal areas and avoided the horizontal discontinuities of the physical processes.Although this approach reduced the horizontal resolution of the physics tendencies, spectral analyses of the kinetic energy of the simulated fields showed that the effective resolution of the dynamic fields was not degraded.
Previous studies have revealed numerous advantages of superparameterization.However, its computational cost has limited the scope of application.Even though the costs of SP-GCMs are lower than high-resolution GCRMs, it is still considerably higher than traditional GCMs, making SP-GCMs less suitable for longer simulations.This problem worsens under a higher-resolution GCM component because the number of CRM domains increases proportionally to the number of GCM columns.Hence, efforts were made to further boost the computational efficiency of SP-GCMs.Mean-state acceleration, proposed by Jones et al. (2015), is a method to significantly speed up SP-GCMs.In mean-state-accelerated models, the time progress of CRM domains is artificially accelerated by a fixed factor, between 2 and 16.If the acceleration factor is set to 2, for example, 1 hr of CRM simulations are coupled to 2 hr of GCM integration.Tendency exchange between the CRM and GCM components is adjusted so that both components run at different time scales consistently.This technique makes possible significantly fast SP-GCM simulations at the cost of artificially shortening the temporal cycles of the solar radiation.Hannah et al. (2020) further boosted their E3SM-MMF model by porting the CRM component to GPUs, taking advantage of the highly parallel nature of the superparameterization.This method is highly promising, considering the overall trends of supercomputing toward GPUs rather than faster CPUs.Although we are currently porting our codes to GPUs as well and preliminary test results are reasonable as expected, details of those computational aspects will be reported in the future.
The goal of this study is to propose a solution to the problem of the computational cost in superparameterization by implementing "blockwise coupling," in which multiple GCM columns are bundled and coupled to one CRM column.Similar modification has been applied to the standard E3SM by Hannah et al. (2021) for the conventional physics parameterization.Similarly, the computational cost of SP-GCMs is expected to be drastically reduced by blockwise coupling because most of the CPU time is consumed by CRM components in SP-GCMs.Furthermore, this modification allows the GCM resolution to be refined without changing the number and size of CRM domains, minimizing the additional computational cost.Because some aspects of the atmosphere (e.g., topographic effects, surface processes, and mesoscale dynamics) are handled by the GCM component in SP-GCMs, increasing the GCM resolution could improve the representation of phenomena sensitive to such GCM-handled processes.In addition, the technique of blockwise coupling will serve as a new tool for the sensitivity experiments of multiscale interactions involving convection.For example, bundling GCM columns to 1,000 × 1,000 km 2 areas for the CRM coupling can disable convective forcing in short wavelengths without affecting the GCM dynamics.The blockwise coupling reduces the horizontal resolution of SP-handled processes including microphysics and radiation, but it is unlikely to severely degrade the overall model resolution because the effective resolution of the dynamics is lower than the nominal grid spacing.Klaver et al. (2020) demonstrated that the effective resolution of CMIP6 spectral models is typically 3-4 times lower than the nominal wavenumber of the spectral truncation.Thus, 4 × 4 physics bundling can possibly help achieve the best match of the effective resolution to that of the dynamics.
In this paper, we introduce SP-MIROC, an SP version of the MIROC6 AGCM (Tatebe et al., 2019) with the aforementioned multiple-to-one blockwise coupling.Model and experimental configurations are described in Section 2. Simulation results of the SP-MIROC is compared against the results of the original MIROC and climatologies in Section 3. Discussions on the effective resolution of blockwise-coupled SP-MIROC is presented in Section 4, and a summary of the work is given in Section 5.
In our companion paper Yamazaki and Miura (2024b), the mechanism responsible for intensified EKWs in SP-GCMs are explored using linear stability analyses with the basic states taken from the MIROC and SP-MIROC experiments described in this paper.

Method
An SP version of a climate model MIROC 6.0 (Tatebe et al., 2019) was implemented by embedding and coupling SCALE-RM (Nishizawa et al., 2015;Nishizawa, Tomita, & SCALE, 2020;Sato et al., 2015), a CRM developed by the RIKEN research institute.Details of the implementation and verification steps of SP-MIROC are given below.

Model
The SP-MIROC consists of MIROC 6.0 AGCM, which serves as the GCM part, and SCALE-RM v5.3.6 as the CRM component.
The dynamic processes in MIROC6, that is, the GCM component, are governed by hydrostatic primitive equations under the T85 horizontal spectral truncation.For physical processes and data I/O, the model has an alternate gridded representation with 256 uniformly spaced zonal grids and 128 Gaussian meridional grids.For coupling to CRM domains, the SP-MIROC uses the same 256 × 128 Gaussian grid.Each grid volume is coupled to one CRM domain in conventional settings of superparameterization, while multiple points from the 256 × 128 Gaussian grid are bundled and coupled to one CRM domain in the blockwise coupling described later.The vertical coordinate in the MIROC6 is the terrain-following σ p hybrid coordinate.There are 40 vertical levels, and the model top is around 7 hPa.In this study, the 81-level configuration of the standard MIROC6 (Tatebe et al., 2019) was not adopted because the radiation scheme used in this study (Sekiguchi & Nakajima, 2008) is designed only for lower altitudes (below 1 hPa), while much higher levels are present in the 81-level configuration.SCALE-RM is a fully compressible nonhydrostatic model and its dynamical core was modified in this study to perform two-dimensional operations instead of three-dimensional ones.The surface elevation is constant within each CRM domain and is averaged from the model topography of the corresponding GCM grids.A fixed zcoordinate is used in SCALE-RM, while the MIROC6 is based on the σ p hybrid coordinate.In this work, the vertical levels of the CRM domains are defined as altitudes of the corresponding GCM levels at the initial time of each SP-MIROC simulation.During model simulations, linear vertical interpolation is performed for GCM-CRM interactions.CRM domains have 32 horizontal grids with 2 km spacing and a periodic condition for the horizontal boundaries.Although the 64-km width of CRM domains is considerably smaller than that of the corresponding GCM columns, the width of CRM domains does not have to match that of the GCM column to reasonably simulate MJOs and can be as narrow as 32 km, as demonstrated by Pritchard et al. (2014).In the current implementation of SP-MIROC, two-dimensional CRM domains are oriented in the zonal direction.Various setups were used in previous studies on CRM orientation: zonal (Grabowski, 2004;Khairoutdinov & Randall, 2001), meridional (Arnold & Randall, 2015), mean wind direction in the lower troposphere (Grabowski, 2004;Tulich, 2015), and a random walk (Tulich & Kiladis, 2021).However, Tulich and Kiladis (2021) found that very different handling of the CRM orientation (4-km mean wind direction and the random walk) resulted in quite similar results.Preliminary experiments in this study also yielded almost identical results from multiple orientation choices.Therefore, the CRM orientation fixed to the zonal direction is unlikely to affect the results of this study.
Most physical parameterizations were deployed only in the CRM domains, in line with the existing SP models.A six-class double-moment bulk microphysics scheme by Seiki and Nakajima (2014) was used despite its relatively high computational cost because a single-moment scheme by Tomita (2008) resulted in significant cold bias in the tropical troposphere.However, the double-moment scheme used in this study (Seiki & Nakajima, 2014) still appears to result in noticeable cloud biases as shown later.The mstrnX radiation scheme (Sekiguchi & Nakajima, 2008) and a Smagorinsky-Lilly type turbulence scheme (Lilly, 1962;Smagorinsky, 1963) were also deployed only in the CRM domains.Some parameterizations were executed exclusively in the GCM component in SP-MIROC.For numerical stability, as described later in this section, the MYNN Level 2.5 planetary boundary layer (PBL) scheme (Mellor & Yamada, 1982;Nakanishi & Niino, 2004) was executed in the GCM component.For simplicity, the surface flux diagnosis and all surface components (e.g., soil, ocean) were processed in the GCM component as in standard MIROC6 AGCM simulations, because they need to be closely coupled to the PBL process.Sea surface temperatures (SSTs) were prescribed, as described in Section 2.3.Land characteristics were prescribed as in standard MIROC6.See Tatebe et al. (2019) for details on data sets used for the land characteristics.
For comparison against SP-MIROC, a simulation using the original version of the MIROC6 AGCM was also performed.This standard MIROC6 experiment will be denoted as CTL-MIROC in this paper.A cumulus parameterization scheme by Chikira and Sugiyama (2010) is deployed in CTL-MIROC.Sea surface temperatures and land characteristics were prescribed identically as in SP-MIROC.See Tatebe et al. (2019) for details of the MIROC6 model, such as other parameterizations.

10.1029/2023MS003836
YAMAZAKI AND MIURA To avoid numerical instabilities, procedures for the time integration in the dynamics of MIROC6, that is, the GCM component of SP-MIROC, was changed.The original MIROC6 employs the leap-frog scheme, which has fair numerical accuracy and low computational cost, for time integration.However, coupling with CRM domains drastically destabilized the leap-frog integration in preliminary experiments because of the additional degrees of freedom.Hence, we adopted the four-stage fourth-order Runge-Kutta method, which is also used in SCALE-RM.Although this increases the computational cost of the GCM dynamics, the cost is negligible compared to that of the CRM domains.

Blockwise CGM-CRM Coupling
In SP-MIROC, one CRM domain is coupled to a horizontal average of 4 × 4 GCM columns, reducing the number of CRM domains to one sixteenth of that in the conventional setup of SP-GCMs.This treatment inevitably reduces the horizontal resolution of physics, including microphysical and radiative heating, but retains the largescale dynamics in the GCM with the original T85 resolution.The resolution of the original physics grid and that of the 4 × 4-bundled grid are compared in Figure 1.While the topography in MIROC is in the original resolution, the 4 × 4-averaged topography is used as the surface altitude in CRM domains.Though most of the coupling procedures were able to operate without issues under this blockwise coupling, redistribution of the surface friction to the atmosphere suffered from a serious problem.If surface winds have a rotating component within 4 × 4-averaging grids of the GCM, averaging the surface stress to transfer to the CRM would cancel out the counter-rotating components of the friction.This allows rotational modes to grow excessively and eventually destabilize the system.To avoid this problem, the surface flux diagnosis and the MYNN Level 2.5 PBL scheme (Mellor & Yamada, 1982;Nakanishi & Niino, 2004) were executed in the GCM component, rather than in the CRM, without the 4 × 4 averaging.
These modifications give rise to constraints in the representation of surface-related properties and processes within the SP-MIROC framework.Due to the flat nature of the topography in the CRM domains, dynamic influences of topography on CRM convection, such as orographic forced asent, cannot be accounted for.On the other hand, large-scale circulations addressed in the GCM component are subject to the influence of GCM topography as in standard GCMs.Another limitation of SP-MIROC pertains to the absence of direct interactions between surface fluxes and the internal circulations of the CRM.While surface fluxes are computed consistently with atmospheric conditions within each GCM column, circulations localized within CRM domains do not exert direct influences on these fluxes.However, it is unlikely that this limitation will severely affect model behaviors because many conventional GCMs also do not incorporate subgrid circulations related to cumulus convection in diagnostics of surface fluxes.Nevertheless, it is desirable to facilitate interactive representation of surface fluxes within the CRM domains in future developments.
In this study, the blockwise coupling is implemented in SP-MIROC but not in standard MIROC6 due to technical difficulties.It is an interesting future topic to implement the blockwise coupling of MIROC6 physics and MIROC6 dynamics, and compare the impact of the blockwise coupling between pure MIROC6 and SP-MIROC.

Coupling Formula
Four variables are coupled between the GCM and CRM components: temperature T, zonal wind u, meridional wind v, and total water mixing ratio q tot , which includes water vapor and all hydrometeors.Labeling those variables in the GCM as Q and those in the CRM as q, tendencies of the superparameterization (∂/∂t) SP are calculated at every timestep of the GCM as shown below: where q indicates the horizontal mean of q within each CRM domain, Q is the horizontal mean of Q within each GCM 4 × 4 column, Δt GCM = 600 s is the GCM time step; (∂/∂t) non SP is the sum of all non-SP tendencies, including dynamics and physics, averaged from the previous GCM timestep to the present.The tendency for the GCM domain (∂Q/∂t) SP is designed to reproduce the internal tendency in the CRM domain (e.g., cloud and radiative forcing) to the GCM, while the tendency for the CRM domain (∂q/∂t) SP forces the horizontal means in the CRM to match those in the GCM columns.To summarize, total tendencies of prognostic variables (∂/∂t) total are written as below: See Appendix A for further details of the coupling formulas in light of discussions by Grabowski (2004).
MIROC6 and SCALE-RM adopt different vertical coordinate systems: σ-p hybrid and fixed height, respectively.Since the altitude of the σ-p layer fluctuates depending on the surface pressure and virtual temperature profile, it is not possible to always match the vertical layers between the GCM and the CRM.Therefore, vertical interpolation is required for all operations between Q and q.To avoid unreasonable overshooting and stencil problems near the boundaries, linear interpolation was adopted in this study.
The linear interpolation, however, has a systematic bias of losing high-order features.The formulation 1 of the SP tendency in the GCM can potentially accumulate such biases, damping out vertical details and disrupting the conservation of variables.To mitigate this problem, the mass-weighted vertical integral of (∂q/∂t) non-SP is interpolated, instead of using the raw values of (∂q/∂t) non-SP .The interpolated integrals are then differentiated on GCM half levels, yielding the interpolated SP-tendency in the GCM columns.Thus, interpolation errors of the integral at a specific altitude only affects the differentiated SP-tendency in the neighboring two layers in opposite signs, and the entire-atmosphere sum is not affected.Nonetheless, the cumulative interpolation scheme presented here is not guaranteed to strictly meet conservation laws of column-integrated water and energy.As noted by Grabowski (2016), a more sophisticated method is desirable to guarantee conservation of variables.

Microphysical Adjustments
Currently, a microphysical adjustment is applied to SP-MIROC to avoid an excessively large bias.In preliminary experiments, SP-MIROC produced large amount of cloud ice in the upper troposphere, significantly reducing outgoing longwave radiation (OLR) and increasing outgoing shortwave radiation (OSR) (figure not shown).
Swapping the microphysics for the single-moment scheme by Tomita (2008) eliminates this bias, indicating the bias is likely caused by the microphysics by Seiki and Nakajima (2014), which has been tuned for regional shortterm simulations, rather than for long-term ones.A modification of the terminal velocity of cloud ice is introduced to partially mitigate the bias.The conversion formula is , where v is the original terminal velocity used by Seiki and Nakajima (2014), and v mod is the modified terminal velocity used in SP-MIROC.This modification is an ad hoc solution to promote the removal of small ice particles in upper-level clouds without affecting larger particles in thick clouds.However, noticeable biases remained in the adjusted SP-MIROC, as shown in Section 3.1.It is an important task for future studies to adjust clouds in a more reasonable manner and eliminate the radiative biases.

Experimental Setup
Five-year simulations were performed using CTL-MIROC and SP-MIROC, and the simulation results were compared from various perspectives.This study focused on simulations with realistic boundary conditions for demonstrating of reproducibility of the blockwise-coupled SP-MIROC.
The SST and sea ice fraction in both CTL-MIROC and SP-MIROC runs were interpolated from the monthly climatology by OISST (Reynolds et al., 2002) from 1981 to 2010.Concentration of greenhouse gases observed from 1995 to 1999 are prescribed, with 1995 as the initial simulation year.For technical simplicity, the initial conditions of all experiments in this study were generated by spinning up a constant-temperature (0°C), motionless atmosphere for 5 years in CTL-MIROC, with the climatological SSTs described above.In the spinup run, variables quickly converged to their respective climatologies in a few months.In SP-MIROC experiments, the GCM component was initialized using the five-year-spinup state introduced above, while the CRM domains were initialized with the values vertically interpolated from corresponding GCM blocks.Because the initial states of the CRM domains described above were horizontally uniform, random noises with an amplitude of 0.01 K was applied to potential temperatures to seed convections.
Furthermore, to check for impacts of the blockwise coupling on SP-MIROC behaviors, additional experiments were performed.Here, the size of GCM blocks for the CRM coupling was changed from 4 to 1 (i.e., no bundling at all), 2, and 8 for sensitivity experiments named 1 × 1-SP, 2 × 2-SP, and 8 × 8-SP, respectively.The simulation period was shortened to 1 year considering the high computational costs.All model configurations other than the block size and integration period remained identical to the standard, 4 × 4-bundled SP-MIROC described earlier, which will be also referred to as 4 × 4-SP.

Comparison Against Observations and Reanalyses
To evaluate the reproducibility of mean states and variability, SP-MIROC and CTL-MIROC results were compared against observations and reanalyses.
To check for the mean-state biases, long-term means were gathered from various sources.The ERA5 reanalysis (Hersbach et al., 2020) was used for wind, humidity, and temperature averaged for 1981-2010.The CERES Energy Balanced and Filled climatologies (Loeb et al., 2018;Wielicki et al., 1996), averaged from 2005 to 2015, were used for the OLR and the OSR.The GPCP precipitation analysis (Adler et al., 2003) was used for mean precipitation rate climatology in 1981-2010.

Wave Composition
To compare the structures of equatorial waves in SP-MIROC and CTL-MIROC and in observations, composite analyses were performed for the EKW, equatorial Rossby wave, mixed Rossby-gravity wave, and n = 1 eastward inertia-gravity wave (EIG) via a procedure similar to that used by Nakamura and Takayabu (2022).Spectral  Wheeler and Kiladis (1999) were used for waves other than the EKW.As for observational composites, the NOAA OLR, the GPCP precipitation, and ERA5 temperatures and winds were used.Nakamura and Takayabu (2022) composited raw variables to the wave phase extracted from OLR; in contrast, the composited variables in this study were prefiltered with the same window as that for OLR to boost the signal-to-noise ratios within the 5-year data period.

Mean States and Variability
Tropospheric zonal mean temperatures, zonal winds and relative humidity in the ERA5 analysis, CTL-MIROC, and SP-MIROC are compared in Figure 2. Overall, the zonal-mean patterns of biases of the two models are similar.CTL-MIROC and SP-MIROC share a warm bias in the polar upper and lower troposphere and a cold bias in lower latitudes.Both models overestimate the core jet speed in the southern troposphere and underestimate its southern extension.Further, both models tend to overestimate the relative humidity throughout the troposphere.In SP-MIROC, the wet bias is larger in the middle troposphere, especially in the tropics.This suggests that the  Figure 3 presents the horizontal distributions of mean precipitation.Figure 3d shows a zonal band of excessive precipitation in the tropical southeast Pacific in CTL-MIROC known as the double-ITCZ.In contrast, the wet bias in the southeast Pacific is reduced in SP-MIROC (Figure 3f).This mitigation of double-ITCZ has been reported in E3SM-MMF (Hannah et al., 2020) as well.Spatial bias patterns of CTL-MIROC and SP-MIROC around the tropical Indian Ocean and tropical western Pacific are similar.In contrast, wet bias of the ITCZ near the Americas in the CTL-MIROC is inverted to a dry bias in SP-MIROC.The global mean precipitation is smaller and closer to the GPCP climatology in SP-MIROC than in CTL-MIROC (Figures 3a,3c,and 3e).However, longwave radiation bias may be responsible for this apparent improvement as discussed later.
The frequency distributions of daily precipitation (Figure 4) reveal that CTL-MIROC produces weak (1-10 mm day 1 ) precipitation too frequently, while this bias is mitigated in SP-MIROC.This feature is common to the improvement in E3SM-MMF (Hannah et al., 2020).Here, because the spatial resolution affects the frequency distributions, 1°data sets (TRMM and GPCP) were converted to 6°fields by bundling and averaging 6 × 6 grid points, and 1.5°outputs by CTL-MIROC were downsampled to 6°fields by averaging 4 × 4-sized blocks.Hence, Figure 4 presents frequency distribution of precipitation averaged in (6°) 2 areas.Although the spatio-temporal mean values differ among data sets (Figure 3), the frequency distributions remain almost unchanged when the individual values are normalized by their respective means (figure not shown).As for very weak precipitation below 0.1 mm day 1 , the frequency is smaller in both CTL-and SP-MIROC than in the TRMM and GPCP estimations.This apparent negative bias is in contrast to E3SM and its SP counterpart, in which such weak precipitation is more frequent than in TRMM and GPCP.However, the sensitivity of the precipitation radar onboard the TRMM satellite is lower for weaker precipitation (Berg et al., 2010), restricting the reliability of the TRMM estimation for faint precipitation especially below 1 mm day 1 .More intensive investigations using observations from highfrequency radar such as the CloudSat Cloud Profiling Radar (Stephens et al., 2008) is desirable for verifying the frequency of such faint precipitation.
Figures 5a and 5b shows that cold cloud ice is globally too abundant in SP-MIROC, while warm water clouds are excessive in CTL-MIROC.Thus, OSR is excessive in both CTL-MIROC and SP-MIROC, while the negative OLR bias is notably larger in SP-MIROC (Figures 5c and 5d).
The abovementioned cloud biases cause insufficient longwave cooling of the atmosphere in SP-MIROC, impacting the global mean atmospheric heat budget:  where SW, LW, LH, SH, TEND, and RES are globally averaged shortwave radiation heating, longwave heating, latent heat release from water, surface sensible heat flux, atmospheric temperature tendency, and residuals, respectively.Values of SW and LW were calculated as net downward fluxes at the top of atmosphere added to the net upward fluxes on the surface.In SP-MIROC, the net longwave cooling is weakened by approximately 20 W m 2 , which is compensated by a reduction of the latent heating proportional to precipitation (Figure 5e).Thus, the apparent mitigation of the excessive global mean precipitation in SP-MIROC (Figure 3) may be a symptom of the cloud bias.

Equatorial Waves
Figure 6 shows the zonal-temporal OLR spectra obtained by the CERES observation, CTL-MIROC, and SP-MIROC in the same manner as that by Wheeler and Kiladis (1999).SP-MIROC excites stronger Kelvin waves than CTL-MIROC, and the former yields an overall spectrum closer to the observation.The intensification of Kelvin waves compared to the standalone GCM in its SP version has been reported by Benedict and Randall (2011) and Hannah et al. (2020) as well.Although antisymmetrical modes are faint in both MIROC models, SP-MIROC produces slightly stronger MRGs and EIGs.Amplitudes of the MJOs are underestimated in both CTL-MIROC and SP-MIROC.Improving representations of the atmosphere-ocean interactions may help reproduce more realistic MJO amplitudes (Benedict & Randall, 2011).
Composite structures of EKWs are shown in Figures 7a-7c.In both CTL-and SP-MIROC, winds and precipitation anomalies are correlated with the OLR phases in a manner similar to that in the analytic data sets.However, the amplitude of CTL-MIROC waves is noticeably smaller for both wind and precipitation perturbations.In contrast, the amplitude of SP-MIROC waves is stronger and closer to that of the analyses.As discussed in Yamazaki and Miura (2024b), top-heavy heating in SP-MIROC, compared to CTL-MIROC, is responsible for the intensified EKWs.Organized precipitating systems in the tropics consist of convective cores and stratiform areas, which have a top-heavy heating profile favorable for EKW growth (Fuchs et al., 2012).Explicit handling of precipitating systems in SP can help reproduce some of the stratiform components, allowing for stronger EKWs.Composites for other equatorial waves are shown in Figures 7d-7l.In CTL-MIROC, the wind amplitude of the equatorial Rossby wave is excessive, while the precipitation anomaly is in a good agreement with the GPCP (Figures 7d-7f).In SP-MIROC, however, the wind amplitude is closer to ERA5, but the precipitation amplitude is too small.CTL-MIROC performs well for MRGs on both wind and precipitation, while precipitation anomalies are too small in SP-MIROC (Figures 7g-7i).With regard to the EIG, SP-MIROC outperforms CTL-MIROC in amplitudes of both winds and precipitation (Figures 7j-7l).These results suggest that the SP-MIROC tends to underestimate the amplitudes of rotational equatorial waves.Sensitivity experiments on bundling block size (see Section 3.3) indicate that the 4 × 4 bundling partially contributes to this bias, while divergent waves (EKWs and EIGs) are unaffected.This suggests that both the blockwise technique and the SP itself may be responsible for the weaker rotational waves.Further investigation, preferably in idealized domains such as aquaplanets for simplicity, will be necessary for reaching detailed conclusions.

Comparison Between Blockwise-Coupled and Conventionally-Coupled SP-MIROC
The implementation of the blockwise GCM-CRM coupling is the most unique feature of SP-MIROC in this study.To validate this treatment, we compared the results of 1 × 1-SP (without the blockwise coupling) and the 4 × 4-SP (with 4 × 4 blockwise coupling).
The mean states obtained from the 1 × 1-SP and the 4 × 4-SP experiments are highly similar for zonal winds, temperatures, relative humidity (Figure 8), mean precipitation (Figure 9), cloud abundance, and top-ofatmosphere radiation (Figure 10).Although SP-GCMs can be sensitive to spatial heterogeneity and produce spurious patterns (Hannah et al., 2020), no unnatural patterns related to the 4 × 4-sized blocks was identified in daily (Figure 11,Figures S1 and S2 in Supporting Information S1) and monthly (Figure S3 in Supporting Information S1) fields.Furthermore, Figure 11 indicates that the simulated synoptic features are highly similar between 1 × 1-SP and 4 × 4-SP, even though 5 days have passed from the common initial condition.Overall, comparison of mean states suggests that the 4 × 4 blockwise coupling does not severely affect the model climatology.
Larger differences between the 1 × 1-SP and the 4 × 4-SP experiments were identified in equatorial waves (Figures 12 and 13).While EKWs are stronger and closer to the observation in 4 × 4-SP, MRGs are stronger in the 1 × 1-SP.Further, although the 4 × 4 bundling is expected to restrict the reproducibility of small features (especially ones with zonal wavenumbers above 20), frequencies of such waves are mostly beyond the datasampling frequency of once per day in this study.Composites of these waves reveal that equatorial Rossby waves and MRGs are better reproduced in the 1 × 1-SP, while EKWs and EIGs are weaker in the 1 × 1-SP.As noted earlier, this finding suggests the blockwise coupling may be slightly unfavorable for rotational waves.However, given the short data period of 1 year, comparisons here contain large statistical uncertainties.Thus, further experiments are needed to reach a decisive conclusion.

Discussions
The horizontal 4 × 4 bundling of GCM columns in SP-MIROC did not result in spurious patterns or severe loss of atmospheric reproducibility (Section 2.1.1).However, there remains the concern of loss of effective resolution because the physics tendencies are smoothed to 4 × 4 GCM grids.Here, the zonal power spectra of kinetic energy, were calculated from daily mean fields.In this study, the analysis exclusively targets the tropics (10°S-10°N) because the zonal size of columns coupled to physics (CRM) in MIROC6 (SP-MIROC) depends on latitudes.While a 4 × 4-bundled grid is roughly 600 km wide in the tropics, it shrinks proportionally to the cosine of the latitude, complicating discussions of effective resolution in higher latitudes.Zonal power spectra were derived for each day and latitude using Fourier transform.Then, the spectra were averaged meridionally between 10°S and 10°N, and temporally for the first year of the simulation, which is the whole period of the sensitivity experiment for block sizes.All MIROC-based outputs were analyzed in the 256 × 128 Gaussian latitude-longitude grid, regardless of the block size, because the GCM component can represent features smaller than the size of blocks coupled to CRM domains.
The kinetic energy in the high wavenumbers increases as the bundled GCM blocks increase (Figure 14), while the small-scale power is noticeably underestimated in the 1 × 1-SP.Although the absolute values of the highwavenumber power differ between the CTL-MIROC and the 4 × 4-SP, no appreciable differences in roll-off wavenumbers can be identified.This suggests that the blockwise technique does not degrade the effective resolution of dynamic variables in SP-MIROC; this result is consistent with the results from the non-SP E3SM (Hannah et al., 2021).As noted in the Introduction, bundling 4 × 4 columns for physics can result in the best match with the dynamics, whose effective resolution is typically 3-4 times lower than the nominal wavenumber of the spectral truncation.This consideration is supported by the fact that the spectra of the 4 × 4-SP typically matches closer to the ERA5 spectra than the other bundling configurations (Figure 14), though the reliability of ERA5 spectra is unclear.These results suggest that the 4 × 4-SP may perform better than the nonbundled SP-MIROC in terms of spectra of GCM-handled variables.However, detailed verification remains a task for future studies.
For variables calculated exclusively by the CRM, such as precipitation and OLR, losses of the highest-wavenumber signals are inevitable.For example, CTL-MIROC and 1 × 1-SP output precipitation every 150 km in the tropics, while 4 × 4-SP calculates it every 600 km.This raises a trade-off between the  spectra of dynamic variables and the small-scale representation of CRM-derived physical variables.For example, considering the extensive loss of high-wavenumber precipitation amplitudes, 8 × 8 bundling may be inappropriate for simulations focusing on regional climates.
When should the blockwise coupling be performed in SP-GCMs, and how many columns should be bundled?GCM dynamics component damps small-scale forcings beyond their effective resolution, whose cutoff scale is typically three to four times larger than their nominal grid spacing.This suggests that 4 × 4 bundling may be a good option regardless of the GCM resolution, although the bundling would limit the applicability of the SP-GCM to smallest GCM disturbances to some extent.In addition, blockwise coupling will be desirable when the assumption of scale separation between GCM-handled features and cumulus convections becomes invalid under high-resolution GCMs.According to the analyses by Arakawa and Wu (2013), the scale separation between cumulus clouds and environment becomes questionable around 30 km.This suggests that the blockwise coupling should be applied to GCM resolutions finer than 30 km, which has already been reached by G. Lin et al. (2022).Grabowski (2016) proposed a large-eddy simulation (LES) version of superparameterization, in which a highresolution GCM is coupled to LES domains.Discussions here suggests that the blockwise coupling would be useful in such LES-SP as well because the horizontal scale of the LES domain would otherwise be close to that of cumulonimbi under GCM resolutions expected in global LES-SP.
As discussed above, while 4 × 4 bundling appeared to perform better in the present study, it is likely influenced by multiple factors such as the effective resolution matching, dissipative properties of the GCM numerics, validity of the scale separation, and other unknown factors, making it difficult to separate the contributions of each aspect.
Further sensitivity experiments with varying GCM resolutions will be needed to distinguish these factors.If the hypotheses made above are to be verified in the future, it will be a good strategy to bundle at least 4 × 4 grids for effective resolution matching; more grids should be bundled so that the resulting blocks are larger than the minimum environmental scale characteristic of their locations, ranging from tens to hundreds of kilometers.

Summary
In this study, the potential scope of superparameterization was extended to higher GCM resolutions by devising and verifying blockwise coupling, in which a horizontal average of multiple GCM columns, instead of one, is coupled to a CRM domain.This modification inevitably reduces the spatial resolution of CRM-handled features in the GCM and limits the applicability of SP-GCMs to small-scale convectively coupled disturbances.However, it drastically reduces the computational cost of SP-GCMs without impacting the overall effective resolution of the model.By coupling SP-GCMs blockwise with higher-resolution GCM components, one can benefit from the finer GCM dynamics without suffering a radical increase of the CRM computational cost.
A blockwise-coupled SP-GCM called SP-MIROC was implemented by coupling the climate model MIROC6 to the CRM SCALE-RM v5.3.6.The blockwise coupling disabled the surface friction for small-scale rotating components of surface winds by spatial averaging, leading to a numerical instability.To solve this problem, the PBL scheme was executed in the GCM without spatial bundling, in contrast to the conventional SP-GCMs.
Tendency transfer from the CRM to the GCM was implemented such that the biases due to vertical interpolation were minimized.Furthermore, the time integration scheme in the MIROC dynamics was replaced by the fourstage fourth-order Runge-Kutta method for addressing the numerical stability.
The implemented 4 × 4-bundled SP-MIROC successfully reproduced horizontal patterns and frequency distributions of precipitation and realistic amplitudes of EKWs, which were underestimated in CTL-MIROC.However, we identified a bias of excessive ice clouds.As discussed in Yamazaki and Miura (2024b) of this study, the amplitude boost of EKWs was enabled by a top-heavy heating in SP-MIROC.No notable artifact was identified in the fields simulated by the 4 × 4-bundled SP-MIROC and the mean states were found to be similar between nonbundled and 4 × 4-bundled SP-MIROC.Comparison of power spectra between the 4 × 4-bundled SP-MIROC and nonbundled SP-MIROC indicated that the effective resolution of dynamic variables was not degraded by the blockwise technique.Rather, spectra in the 4 × 4-bundled SP-MIROC were more realistic than those in the nonbundled SP-MIROC.
This study leads to multiple interesting topics and issues left for future works.One topic arises from the fact that the power spectra of 4 × 4-bundled SP-MIROC is more realistic than the nonbundled one.This raises a hypothesis that the CRM domains of SP "should" be coupled blockwise, possibly at the effective resolution of version of the climate model MIROC6 named superparameterized (SP)-MIROC was developed • Computational cost is reduced by coupling a cloud-resolving model domain to a block of multiple global climate model columns, instead of one • Blockwise-coupled SP-MIROC reproduced most benefits of SP without losing the effective resolution of dynamic variables Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.Land fraction and topography in the (a, c) the original physics grid in the T85 MIROC6 and (b, d) the 4 × 4-bundled grid used for cloud-resolving model coupling in this study.

Figure 2 .
Figure 2. Upper panels: zonal mean (a) temperature, (b) zonal wind, and (c) relative humidity of the ERA5 climatology.Middle panels: zonal and 5-year mean (d) temperature, (e) zonal wind, and (f) relative humidity in CTL-MIROC represented by contours and their differences from ERA5 climatologies (shading).Bottom panels: same as the middle panels, but for SP-MIROC.Differences are significant at the 95% confidence level where stippled.

Figure 3 .
Figure 3. Upper panels: (a) annual mean precipitation by the GPCP climatology and (b) its zonal-mean values (black line) and the corresponding values simulated by CTL-MIROC (blue line) and SP-MIROC (red line).Strippling indicates biases significant at the 95% confidence level.Middle panels: (c) annual mean precipitation in CTL-MIROC and (d) its difference from the GPCP climatology.Lower panels: same as the middle panels, but for SP-MIROC.While left panels (a, c, and e) are drawn in the original resolution of each data set, precipitation biases (d, f) are computed and drawn commonly in the 4 × 4-bundled MIROC grid.

Figure 4 .
Figure 4. Frequency distribution of daily precipitation between 50°S and 50°N.Precipitation fields were downsampled to 6°grids before deriving the frequency distribution.

Figure 6 .
Figure 6.Upper panels: spatio-temporal spectra of the equator-symmetrical components of daily outgoing longwave radiation (OLR) by the (a) NOAA OLR observation, (b) CTL-MIROC, and (c) SP-MIROC.Five-year fields between 15°S and 15°N are used.Gray lines indicate dispersion relations of equatorial Kelvin waves, equatorial Rossby waves, and symmetric inertia-gravity waves at equivalent depths of 12, 25, and 50 m.Lower panels: the same as the upper panels but for antisymmetrical components and dispersion relations of MRGs and antisymmetric inertia-gravity waves.

Figure 14 .
Figure 14.Zonal spectra of kinetic energy at 200 hPa in ERA5, CTL-MIROC, and SP-MIROC experiments.In SP-MIROC experiments, the spectra were derived from the large-scale winds handled by the global climate model (GCM) component.The spectra are averaged meridionally for 10°S-10°N and temporally for the first year of the simulations.The Nyquist wavenumbers corresponding to sizes of CRM-coupled blocks are 128 (1 × 1-SP), 64 (2 × 2-SP), 32 (4 × 4-SP), and 16 (8 × 8-SP), respectively.However, the GCM dynamics handles large-scale winds up to a zonal wavenumber of 85, regardless of the block sizes in SP-MIROC.

Modeling Earth Systems
this study are listed in Table 1.To avoid complexities arising from potentially wavenumber-dependent excitation mechanisms, a narrow window employed by Nakamura and Takayabu (2022) was adopted for the EKW; the mechanisms are explored in details in Yamazaki and Miura (2024b) of our work.Because the excitation mechanisms of other waves are beyond the scope of this work, broadband windows employed by

Table 1
Filters Used for the Composite Analyses of Equatorial Waves

Journal of Advances in Modeling Earth Systems
YAMAZAKI AND MIURA