The Nonhydrostatic ICosahedral Atmospheric Model for CMIP6 HighResMIP simulations (NICAM16-S): experimental design, model description, and impacts of model updates

The Nonhydrostatic ICosahedral Atmospheric Model (NICAM), a global model with an icosahedral grid system, has been under development for nearly two decades. This paper describes NICAM16-S, the latest stable version of NICAM (NICAM.16), modified for the Coupled Model Intercomparison Project Phase 6, High Resolution Model Intercomparison Project (HighResMIP). Major updates of NICAM.12, a previous version used for climate simulations, included updates of the cloud microphysics scheme and land surface model, introduction of natural and anthropogenic aerosols and a subgrid-scale orographic gravity wave drag scheme, and improvement of the coupling between the cloud microphysics and the radiation schemes. External forcings were updated to follow the protocol of the HighResMIP. A series of short-term sensitivity experiments were performed to determine and understand the impacts of these various model updates on the simulated mean states. The NICAM16-S simulations demonstrated improvements in the ice water content, high cloud amount, surface air temperature over the Arctic region, location and strength of zonal mean subtropical jet, and shortwave radiation over Africa and South Asia. Some long-standing biases, such as the double intertropical convergence zone and smaller low cloud amount, still exist or are even worse in some cases, suggesting further necessity for understanding their mechanisms, upgrading schemes and parameter settings, and enhancing horizontal and vertical resolutions.

1 Introduction

Moist processes play a crucial role in the formation of the earth's climate. Moist convection redistributes mass, energy, and momentum of the atmosphere to form large-scale circulation. Clouds are coupled with large-scale circulation through latent and radiative heating, which can affect climate sensitivity. The accurate treatment of such interactions between clouds and circulation requires high-resolution global cloud-resolving models (Bony et al., 2015; Satoh et al., 2019).

This, as well as an increasing demand from society to project extreme weather, such as tropical cyclones, motivated us to perform climate simulations of the present day and future using the Nonhydrostatic ICosahedral Atmospheric Model (NICAM) (Tomita and Satoh, 2004; Satoh et al., 2008, 2014) with a 14 km mesh. Kodama et al. (2015) and Satoh et al. (2015) provided brief descriptions of the model (hereafter referred to as NICAM.12) and experimental design of the climate simulations. This unique dataset of the high-resolution climate simulation, whose overall performance was reported in Kodama et al. (2015), has been used in many studies with a focus on tropical cyclones (Satoh et al., 2015, 2018; Yamada et al., 2017, 2019; Matsuoka et al., 2018; Sugi et al., 2020), extratropical cyclones (Kodama et al., 2019; McCoy et al., 2019; Satoh et al., 2018), intraseasonal oscillations, such as the Madden–Julian Oscillation (MJO) (Kikuchi et al., 2017; Nakano and Kikuchi, 2019), tropical synoptic-scale waves (Fukutomi et al., 2016), cloud radiative feedback (Chen et al., 2016; Satoh et al., 2018; Noda et al., 2019), and regional to global precipitation (Kilpatrick et al., 2017; Satoh et al., 2018; Na et al., 2020; Takahashi et al., 2020).

Also, some significant climate biases have been identified in the simulation (Kodama et al., 2015), and great effort has been devoted to improving the model for better performance of the simulated climate in a physically consistent manner. Major updates from NICAM.12 to NICAM.16, a stable version of NICAM released in 2017, included the cloud microphysics scheme based on a comparison with the Tropical Rainfall Measuring Mission (TRMM) satellite observations (Roh and Satoh, 2014; Roh et al., 2017), introduction of a wetland scheme in the land surface model (Nitta et al., 2017), implementation of the coupling between cloud microphysics and radiation that considers the nonsphericity of ice particles (Seiki et al., 2014), and introduction of a subgrid-scale orographic gravity wave drag scheme. In addition, some parameters related to the processes of surface albedo and sea ice have been revised. These model updates generally reduce the biases of the simulated mean states, as discussed below. NICAM.16 has been further modified to support the external forcings of natural and anthropogenic aerosols and the solar cycle defined in the Coupled Model Intercomparison Project Phase 6 (CMIP6) High Resolution Model Intercomparison Project (HighResMIP) protocol (Haarsma et al., 2016). This special version of NICAM.16 for HighResMIP is labeled NICAM16-S, where “-S” represents the use of a single-moment cloud microphysics scheme. A double-moment cloud microphysics scheme is also available in NICAM.16 (Seiki and Nakajima, 2014; Seiki et al., 2014, 2015b; Satoh et al., 2018). However, the double-moment scheme was not used for the HighResMIP simulations and hence is not described in this paper. NICAM.16 is not coupled with an ocean model, and the DECK (Diagnostic, Evaluation and Characterization of Klima) and CMIP historical simulations (Eyring et al., 2016) are not presented in this study. Note that a coupled ocean–atmosphere model, NICAM-COCO (Miyakawa et al., 2017), is being developed.

This section has provided a summary of NICAM16-S with a focus on its differences from NICAM.12. Section 2 of this paper presents the experimental design of the HighResMIP simulations and a series of sensitivity experiments by NICAM16-S. Section 3 explains the detailed model updates of NICAM16-S from NICAM.12 and their impacts on the simulated mean states. Section 4 briefly describes the resolution dependency of NICAM16-S. Section 5 reports on the computational aspects of the simulation. Section 6 provides a brief summary of the comparison results.

2 Experimental design

2.1 Spatial and temporal resolutions

Three sets of model configurations were prepared for the HighResMIP simulations, and initial and boundary conditions were made for each. NICAM16-S models with specific horizontal resolutions were formally labeled NICAM16-7S (56 km mesh), NICAM16-8S (28 km mesh), and NICAM16-9S (14 km mesh) in CMIP6 as source IDs. The horizontal mesh size is evaluated as a square root of the mean area of each grid cell (Satoh et al., 2014). The number n in NICAM16-nS is a grid division level (glevel), which denotes the number of subdivisions of the icosahedron to generate a mesh (Tomita et al., 2001). The physics schemes, including parameters, and the initial and boundary conditions are common among different horizontal resolutions except for those explicitly noted in Sects. 2 and 3.

The number of vertical levels is 38, with a model top height of around 40 km, equivalent to the previous climate simulations (Kodama et al., 2015). The interval between each vertical layer increases from 160 m to 2 km as the altitude increases from the ground to 25 km (see the “K38” setting in Fig. 1 of Ohno et al., 2019). Even at such a low vertical resolution, atmospheric phenomena of interest, including tropical cyclones, the MJO, and the diurnal precipitation cycle, may be simulated practically and accurately, as we have confirmed in the previous study using a 14 km horizontal mesh (Kodama et al., 2015). As a caveat, a coarse vertical resolution in the upper atmosphere leads to an overestimation of the cirrus cloud amount (Seiki et al., 2015b; Ohno et al., 2019) and may cause a different response of high cloud to warmer climate (Ohno et al., 2019). Also, it has been suggested that the vertical resolution should be increased when the horizontal resolution is increased in terms of atmospheric gravity waves (Lindzen and Fox-Rabinovitz, 1989; Polichtchouk et al., 2019). Such coarse vertical resolution in this study could change vertical propagation of gravity waves and zonal wind in the stratosphere (Watanabe et al., 2015; Skamarock et al., 2019).

The time step of the dynamics (Δt in Satoh et al., 2008) was set to 4, 2, and 1 min in NICAM16-7S, NICAM16-8S, and NICAM16-9S, respectively. The time loop in the model is based on the dynamics, and physics schemes with a time step smaller or greater than that of the dynamics are subcycled or skipped, as appropriate. Specifically, a time step of 30 s was used in the cloud microphysics scheme in NICAM16-7S, NICAM16-8S, and NICAM16-9S. A time step of 1 min was used in the turbulence (mainly for the planetary boundary layer) and land and ocean surface schemes in NICAM16-7S, NICAM16-8S, and NICAM16-9S. The radiation scheme, which requires considerable computational time, was executed every 40, 20, and 10 min in NICAM16-7S, NICAM16-8S, and NICAM16-9S, respectively. The gravity wave drag scheme was called at the same time step as that of the dynamics.

2.2 HighResMIP simulations and sensitivity experiments

Table 1 shows the integration periods for the HighResMIP simulations. For the Tier 1 simulations, which start from 1 January 1950, the initial condition of the atmosphere was taken from the ERA-20C reanalysis (Poli et al., 2016). The NICAM16-7S and NICAM-8S simulations strictly followed the HighResMIP protocol and continued until 31 December 2014. The HighResMIP Tier 3 simulations using NICAM16-7S and NICAM16-8S started on 1 January 2015 as a continuation of the Tier 1 simulations and ended on 31 December 2050. The high computational cost hindered us from running NICAM16-9S for 100 years and thus a time slice approach was adopted instead. Specifically, climate simulations were performed in the following timeframes: 1950–1960, 2000–2010, and 2040–2050. For the simulations starting from 1 January 1950 or 1 January 2000, the initial land conditions prescribed for NICAM16-7S, NICAM16-8S, and NICAM16-9S were taken from the monthly mean climatology of the simulation by NICAM with a mesh size of 220 km (glevel 5) under present-day conditions. It was performed for 10 years, and the last 5 years of data were used to obtain the monthly mean land climatology. This approach is the same as that used in previous climate simulations (Kodama et al., 2015). This could partly reduce the initial shock of the land surface model, even though it may cost more than several years for some land variables, such as soil moisture, to fully settle down (not shown). The initial land condition for the future time slice run with NICAM16-9S was obtained by interpolating the output of the NICAM16-8S simulation.

Table 1 List of HighResMIP simulations.

Table 2 List of sensitivity experiments.

In addition to the formal HighResMIP simulations, we performed a series of short-term sensitivity experiments to evaluate the impacts of the model changes on the simulated climatology, as listed in Table 2. Here, the model configuration of the REFFIX run is equivalent to that used in the formal HighResMIP simulations. As noted in Sect. 3.7, we often prefer to use a slab ocean model with nudging toward the boundary sea surface temperature (SST) rather than the fixed SST condition requested by the HighResMIP protocol because of better performance in the simulated precipitation pattern (Kodama et al., 2015), particularly with a horizontal mesh size of 14 km (Sect. 3.7). Therefore, both the fixed SST and slab ocean configurations (REFFIX and REFSLB runs, respectively) were tested in the sensitivity experiments, and we used the REFFIX run with 56 km mesh and the REFSLB run with 14 km mesh as the reference (REF) runs for the other sensitivity experiments. Impacts of the model updates described in Sect. 3 on the simulated climate states were individually tested by switching off each update. Additionally, the DDT2M, DDT1M, RDT20M, and RDT10M runs were performed to check the sensitivity of the simulated climate to the time step of the model. These sensitivity experiments were started on 1 June 2004 and integrated for 1 year. An exception was the NOLND and the REF runs, which were performed for 4 years to make the land surface state settle down. The initial date was chosen to ensure consistency with previous NICAM studies (e.g., Kodama et al., 2012; Seiki et al., 2015a; Noda et al., 2016). An integration period of 1 year, or even less, is sufficient to evaluate the basic state of the atmosphere, such as cloud amount, precipitation, radiation, and temperature (e.g., Phillips et al., 2004; Noda et al., 2010; Kodama et al., 2012; Williams et al., 2013; Miyakawa et al., 2018; Miyakawa and Miura, 2019; Stevens et al., 2019; Hohenegger et al., 2020), and tropical variability, including diurnal cycle, tropical cyclones, and MJO (e.g., Sato et al., 2009; Kinter et al., 2013; Stevens et al., 2019; Matsugishi et al., 2020). The interannual variability of the HighResMIP simulation by NICAM16-7S (Table 1) over 101 years was diagnosed to distinguish the impacts of the model changes from internal variability in a rough manner. The simulation data were regridded to the same grid of observations (Table 3) unless it is explicitly specified.

Table 3 List of observational datasets.

2.3 External forcings and boundary conditions

External forcings and boundary conditions of the simulations followed the HighResMIP protocol (Haarsma et al., 2016). Historical and SSP5-8.5-scenario settings (O'Neill et al., 2016) were used in the Tier 1 and 3 simulations, respectively.

Daily quarter-degree SST and sea ice mass (ICE) prescribed for the model were obtained from HadISST 2.2.0.0 (Kennedy et al., 2017). The SST dataset was extended from 2016 to 2050 using the trend obtained from a CMIP5 model ensemble mean following the RCP8.5 scenario and historical variability from 1980 to 2015 (Kennedy et al., 2019). Because HadISST provides historical sea ice concentration (SIC), ICE was diagnosed from SIC for NICAM (see Sect. 3.7). The future SIC was estimated from the future SST data and the observed relationship between SST and SIC (HighResMIP, 2020; Kennedy et al., 2019). Both the SST and ICE were fixed to the boundary conditions in the HighResMIP Tier 1 and 3 simulations.

Figures 1 and 2 display the decadal mean SST and ICE, respectively, prescribed for the model; Fig. 3 exhibits their global mean variability. Greater warming over the maritime continent, the Indian Ocean, and the edge of the polar regions are found in the 2000s compared with the 1950s, whereas cooling is noticed in the North Pacific and the North Atlantic. The SST in the 2040s has larger values almost everywhere compared with that in the 2000s, especially in the midlatitudes, the equatorial eastern Pacific, the tropical Atlantic Ocean, and the edge of the Arctic regions. Similar tendencies are observed for the distribution of the SST trend (not shown). The global mean SST is 17.9 ^∘C in the 1950s, 18.2 ^∘C in the 2000s, and 19.0 ^∘C in the 2040s. ICE continues to decrease from the past to the future. The global mean ICE is nearly halved by the 2040s compared with that in the 1950s.

Figure 1 Decadal mean horizontal distribution of sea surface temperature (SST) prescribed in the model averaged for the 1950s (a). Differences between the 2000s and the 1950s (b), the 2040s and the 1950s (c), and the 2040s and the 2000s (d) are also shown. Units are ^∘C.

Figure 2 Same as Fig. 1 but for sea ice mass (ICE) in kg m⁻².

The global annual mean of greenhouse gas (GHG) concentrations (Meinshausen and Vogel, 2016; Meinshausen et al., 2017; Meinshausen and Nicholls, 2018) was prescribed in the model. Specifically, CO₂, CH₄, N₂O, chlorofluorocarbons (CFC-12, CFC-11, CFC-113, CFC-114, and CFC-115), hydrochlorofluorocarbons (HCFC-22, HCFC-141b, and HCFC-142b), hydrofluorocarbons (HFC-134a, HFC-32, HFC-125, HFC-143a, and HFC-152a), CCl₄, CF₄, SF₆, and C₂F₆ were considered as GHGs in the model. In addition, historical and future monthly concentrations of the three-dimensional ozone field (Hegglin et al., 2016, 2018) was prescribed in the model.

Figure 3 Global mean SST (a; in ^∘C) and ICE (b; in kg m⁻²) prescribed in the model. The annual mean (thin curve) and decadal running mean (thick curve) are shown.

Mass and number concentrations of natural aerosols prescribed in the model were obtained from a low-resolution NICAM simulation with an online aerosol module based on the Spectral Radiation-Transport Model for Aerosol Species (SPRINTARS) (Takemura et al., 2000, 2002, 2005, 2009; Goto et al., 2008, 2011). The simulated climatology of aerosols in NICAM was previously validated (Goto et al., 2018; Suzuki et al., 2008). For HighResMIP, NICAM with a mesh size of 220 km was performed for 100 years using natural aerosol emissions with the anthropogenic aerosol module MACv2-SP (Fiedler et al., 2017; Stevens et al., 2017; see Sect. 3.4) under a perpetual 2012 condition, and the data for the last 90 years were averaged to obtain a monthly mean climatology of aerosol mass and number concentrations of cloud condensation nuclei (CCN) from a natural origin. A lower-bound limiter of 50 cm⁻³ was applied to the CCN prescribed for the model to avoid numerical instability in the cloud microphysics scheme. Figure 4a shows the annual mean of the natural aerosol optical thickness simulated by NICAM16-7S with the natural aerosol mass and CCN (Fig. 4b) prescribed in the model.¹ The climatology of natural aerosol mass and CCN is invariant year by year, whereas anthropogenic aerosols from MACv2-SP are time dependent in the historical and future simulations (Fiedler et al., 2019; Stevens et al., 2017). Further, the extinction coefficient, the single-scattering albedo, and the asymmetric factor were overwritten with the stratospheric aerosol dataset (Thomason et al., 2015) above the tropopause to introduce the effect of volcanic eruptions on the radiation field in a consistent way among different models participating in the HighResMIP.

Figure 4 (a) Annual mean natural aerosol optical thickness averaged for June 2004–May 2005 as simulated by NICAM16-7S (NOANTAER run in Table 2). The data are gridded at 0.56^∘ in longitude and latitude. (b) The annual mean number concentration of cloud condensation nuclei (CCN) from the natural origin at 1 km above sea level prescribed for the model. The data are gridded at 2^∘ in longitude and latitude and the unit is cm⁻³. The lower bound of CCN, 50 cm⁻³ (Sect. 2.3), is shown in white (no shading).

Similar to the implementation of MIROC6 (Tatebe et al., 2019), historical monthly mean solar forcings (Matthes et al., 2017b, a) were prescribed as total solar irradiance and solar irradiance spectra in the radiation scheme mstrnX (Sekiguchi and Nakajima, 2008). In terms of land surface processes, a monthly climatology (2004–2013) of leaf area index was obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) product (MCD15A2.005, Shabanov et al., 2005; Yang et al., 2006).

As briefly noted in Satoh et al. (2014), a spatial filter was applied to smooth the model topography to avoid numerical instability. Specifically, a hyper-diffusion was repeatedly applied to the GTOPO30 (https://doi.org/10.5066/F7DF6PQS, last access: 20 January 2021), a global digital elevation model with a horizontal spacing of approximately 1 km, to meet a specific criterion for the maximum elevation gradient. The maximum elevation gradient was set to 0.01, 0.01414, and 0.02 m m⁻¹ for NICAM16-7S, NICAM16-8S, and NICAM16-9S, respectively, and the resulting topography is called “A-topography”. Note that in previous NICAM studies using a 14 km horizontal mesh (e.g., Kodama et al., 2015), “B-topography”, in which A-topography with a 28 km mesh was interpolated to 14 km mesh grid points, was often used for the sake of stable integration. For HighResMIP, A-topography was used to better represent steeper mountains and their effects on atmospheric phenomena.

3 Model description and impact of model updates on the simulated fields

3.1 Overview

The dynamical core and numerical filters in NICAM16-S are the same as those in NICAM.12. NICAM adopts a fully compressible nonhydrostatic system as the governing equations of the dynamics (Tomita and Satoh, 2004; Satoh et al., 2008). The horizontal discretization is an icosahedral grid system modified with spring dynamics for homogeneity on the sphere (Tomita et al., 2002). Divergence damping and second-order Laplacian horizontal diffusion are used to stabilize the integration (Satoh et al., 2008). Additionally, first-order Laplacian horizontal diffusion is applied above an altitude of 20 km to avoid spurious wave reflection at the model top.

Table 4 gives a summary of the physics schemes used in NICAM16-S and NICAM.12. A single-moment bulk cloud microphysics scheme, NICAM Single-Moment Water 6 (NSW6), that solves mass concentrations for six water categories – vapor, cloud water, cloud ice, rain, snow, and graupel (Tomita, 2008; Roh and Satoh, 2014; Roh et al., 2017) – is used instead of a combination of convection and large-scale condensation schemes. While most climate models use convection and large-scale condensation schemes, even for a mesh size around 14 km, we used the cloud microphysics scheme to represent interactions between clouds and circulation in an explicit way. This not only lowers the cost of model development but also reduces the uncertainty in the results arising from highly arbitrary tuning. This approach has also been tested by other researchers besides the NICAM group (Maher et al., 2018; Hohenegger et al., 2020). Global mean precipitation is constrained by radiative cooling in large-scale clear-sky regions, which can be captured by a model with relatively coarse resolution without the convection and large-scale condensation schemes. The simulated climatology of the precipitation pattern, even with the lowest resolution setting (NICAM16-7S), is comparable with the observed patterns, as shown below, although our choice leads to patchy precipitation behavior and dry/wet bias in the middle/lower troposphere in the simulation (Miyakawa et al., 2018). Similar precipitation behavior was also reported in a climate model study with a mesh size of around O(10² km) without a convection scheme (Maher et al., 2018). In terms of clouds, Seiki et al. (2015b) conducted NICAM simulations with 28 and 14 km meshes to study the impact of the vertical resolution on the simulated cirrus clouds. They found a similar vertical resolution dependency between the 14 and 28 km meshes. Ohno et al. (2019) used a 28 km mesh NICAM and found a high cloud response to SST increases that is comparable with results using 7 and 14 km meshes (Iga et al., 2007). In terms of MJO, Takasuka et al. (2018) performed an aqua-planet experiment with a 56 km mesh NICAM to investigate MJO-like disturbances. Yoshizaki et al. (2012) and Takasuka et al. (2015) even performed NICAM with a mesh size larger than 100 km without the convection and large-scale condensation schemes and found MJO-like disturbances in the simulation.

Table 4 Physics schemes in NICAM16-S and NICAM.12.

A modified version of the Mellor–Yamada level 2 scheme (Nakanishi and Niino, 2006; Noda et al., 2010) is used to simulate the planetary boundary layer. The radiation scheme, mstrnX (Sekiguchi and Nakajima, 2008), is a broadband model with 29 radiation bands as used here. The land surface model, Minimal Advanced Treatments of Surface Interaction and RunOff (MATSIRO) (Takata et al., 2003) solves land states such as soil temperature, soil moisture, and land surface fluxes. Ocean surface fluxes are calculated following Louis (1979) with a modified roughness length for strong surface wind conditions (Fairall et al., 2003; Moon et al., 2007). The conventional orographic gravity wave drag scheme (McFarlane, 1987) is used to introduce the effect of vertically propagating subgrid-scale orographic gravity waves on the momentum tendency of the atmosphere. Though we did not fine-tune the model due to heavy computational cost, we crudely turned parameters of sea ice thickness with NICAM16-7S (Sect. 3.7) and the gravity wave drag scheme with NICAM16-9S (Sect. 3.8). We did not return the model at each resolution under the principles of the HighResMIP.

As we discuss in Sect. 2.2 and Table 2, a series of short-term sensitivity experiments were performed to monitor the impacts of several model updates on the simulated climatology. Figure 5 (right part of each panel) summarizes the impacts of the model changes on the global mean climate. All the significant impacts of the model changes shown here could be qualitatively reproduced even when the analysis period was limited to the last 6 months (not shown). The REFFIX and REFSLB runs with each horizontal mesh and the observations are shown on the left in each panel in Fig. 5. We will discuss these impacts along with the details of the model updates later in this section.

Figure 5 Global annual means of surface air temperature (a), precipitation (b), top-of-atmosphere (TOA) outgoing longwave radiation (OLR) (c), TOA outgoing shortwave radiation (OSR) (d), ice water path (e), liquid water path (f), high cloud amount (g), middle cloud amount (h), low cloud amount (i), surface net downward longwave radiation (j), surface net downward shortwave radiation (k), surface latent heat flux (l), and surface sensible heat flux (m). The global annual means are averaged over June 2004–May 2005, and their units are shown at the upper-right corner of each panel. Blue shading shows the interannual variability (2σ, detrended) estimated from the HighResMIP NICAM16-7S run over 1950–2050 (Table 1). The left part of each panel contains plots of global annual means simulated by NICAM16-7S (56 km mesh; blue), NICAM16-8S (28 km mesh; green), and NICAM16-9S (14 km mesh; red), which were performed under the fixed SST condition (filled circles; REFFIX run in Table 2) and with the slab ocean condition (filled squares; REFSLB run in Table 2). Blue and red lines are the reference (REF) runs with 56 km mesh and 14 km mesh, respectively. Observational values taken from the JRA-55 reanalysis (surface air temperature), GPCP (precipitation), CERES (radiation), and ISCCP (cloud amount) are shown as gray lines. The right part of each panel shows differences between the REF run and each sensitivity run (NOCLD, NONS, NONSI NOAER, NOALB, NOSIC, and NOGWD runs in Table 2). simulated by NICAM16-7S (56 km mesh; blue multiplication sign), and NICAM16-9S (14 km mesh; red multiplication sign). Results outside the scale range are shown with digits.

3.2 Cloud microphysics

Recently, Roh and Satoh (2014) and Roh et al. (2017) significantly revised the NSW6 scheme based on comparisons with TRMM observations. We used the revised version of the NSW6 scheme in NICAM16-S and found improvements in the simulated climatology, as shown below.

The NSW6 scheme originated with Lin et al. (1983) and Rutledge and Hobbs (1984). Tomita (2008) modified their cloud microphysics scheme to ensure consistency with the thermodynamics used in NICAM. Tomita (2008) also simplified it to reduce the calculation cost for high-resolution global simulations. NSW6 was evaluated by comparing the simulated optical properties with satellite observations (Satoh et al., 2010; Kodama et al., 2012; Hashino et al., 2013, 2016; Roh and Satoh, 2014, 2018; Roh et al., 2017) using satellite simulators, specifically, the CFMIP Observational Simulator Package (Haynes et al., 2007; Chepfer et al., 2008; Bodas-Salcedo et al., 2008, 2011) and Joint Simulator (Matsui et al., 2009; Masunaga et al., 2010; Hashino et al., 2013). It was revised in each stage of the version management of NICAM (Kodama et al., 2012; Roh and Satoh, 2014; Roh et al., 2017). The revision of the NSW6 scheme by Roh and Satoh (2014) and Roh et al. (2017) represents a significant change between NICAM.12 and NICAM16-S.

Table 5 summarizes key changes in the NSW6 scheme by Roh and Satoh (2014) and Roh et al. (2017). In short, the revision aimed to enhance organizations of tropical convective cloud systems by assuming lighter precipitation of graupel and snow and by moderating the development of cloud ice. Finally, they successfully reproduced the vertical structures of shallow, congestus, and deep convective clouds over the tropics compared to TRMM and CloudSat satellite observations. They also used microwave satellite observations to evaluate the simulated results (Roh and Satoh, 2018). Therefore, improvements in tropical cloud systems with the revised scheme are robust in terms of optical signals (see the original papers for more details). Note that the separation of convective and stratiform systems in Roh and Satoh (2014) was omitted in this study because of its small impact (not shown) despite the high computational cost.

Table 5 Summary of the key changes in the NSW6 scheme by Roh and Satoh (2014) and Roh et al. (2017).^*

^* The particle size distribution for rain was also revised in the original paper, but the revision is not used in the latest version because it had a small impact and reduced the computational efficiency. In addition, the assumption that cloud ice does not precipitate is inconsistent with some of the other ice cloud microphysics assumptions, but it is used to tune the model to the observed high cloud signals over the tropics.

The sensitivity experiments with and without the update of the cloud microphysics scheme (the REF and NOCLD runs, respectively) were compared. The most noticeable impact was an increase in ice water content (IWC) by more than 2-fold (Fig. 5e), and this mostly accounts for the snow category in the cloud microphysics scheme. Figure 6 shows the meridional-height cross section of the observed and simulated zonal mean IWC and its breakdown into the categories of cloud ice, snow, and graupel in the NOCLD and REF runs. The simulated IWC is largely underestimated using the model without the updated cloud microphysics scheme (the NOCLD run), as also shown in Seiki et al. (2015a), and it becomes comparable with CloudSat observations with the update (the REF run). A noticeable increase in the snow category was seen in the tropical upper troposphere and midlatitude storm-track region. Cloud ice is also increased in the upper troposphere. Graupel increased in the tropical middle troposphere but decreased in the storm-track region. As a result, global mean column-integrated cloud ice and snow increased by 24 % and 399 %, respectively, and that of graupel decreased by 8.4 %. The net increase in IWC is consistent with the decelerated development of IWC by the modified mass and diameter relationship of snow, which reduces the snow density (Table 5, d) by the diminished efficiency of accretion of cloud ice by snow (Table 5, g) and by ignoring the accretion of snow and cloud ice by graupel (Table 5, f). The increase/decrease in graupel by the model update may be linked to the differences in cloud processes (convection versus synoptic systems).

Figure 6 Annual mean of the zonal mean ice water content (10⁻⁶ kg m⁻³) in CloudSat observations (a) and the NOCLD (b) and REF (c) runs by NICAM16-9S. The breakdown of the simulated ice water content (IWC) into cloud ice (d), snow (e), and graupel (f) is shown in the bottom panels. The contours show the NOCLD run and shading shows the difference between the REF and NOCLD runs. The analysis data are 0.25^∘ (a) and 2.5^∘ (b–f) gridded data, and the vertical axis is the altitude in km above sea level.

Despite the drastic increase in mass concentrations of snow and cloud ice, the amount of high clouds, particularly, optically thin clouds, is reduced by the update (Fig. 5g). Consistently, global mean outgoing longwave radiation (OLR) is increased by about 4 W m⁻² (Fig. 5c), opposite to that found in Roh et al. (2017). In addition, the decrease in top-of-atmosphere (TOA) brightness temperature by the update is very small (Fig. 7) compared with that in Roh et al. (2017). These differences between Roh et al. (2017) and this study can be mostly explained by the new coupling procedure between cloud microphysics and radiative transfer, as described in Sect. 3.3, and partially explained by the different treatment of the terminal velocity of cloud ice in the reference runs. Eventually, the cloud ice terminal velocity was set to zero in both Roh et al. (2017) and this study. Unlike this study, Roh et al. (2017) performed their reference run with nonzero cloud ice terminal velocity, as diagnosed by Heymsfield and Donner (1990), and their comparison before and after the scheme update includes the effects of reduced cloud ice terminal velocity. The reduction of the cloud ice descent speed, as in Roh et al. (2017), could increase and elevate high clouds and decrease the OLR (Kodama et al., 2012). The low cloud amount increased (Fig. 5i) as a result of the compensation between an increase in medium and thick clouds and a decrease in thin clouds. These results indicate that the clouds grow thicker on average by updating the cloud microphysics scheme in this study.

Figure 7 TOA brightness temperature near 11 µm at 00:00 UTC 6 June 2004 from the GridSat product (a) and the NOCLD (b) and REF (c) runs by NICAM16-9S. The display style follows Fig. 1 of Roh et al. (2017). The grid interval in (b) and (c) is 0.14^∘.

3.3 Coupling between cloud microphysics and radiative transfer

In NICAM16-S, the cloud microphysics scheme is fully coupled with the radiation scheme mstrnX (Sekiguchi and Nakajima, 2008). The effective radius of hydrometeors is calculated with the same assumption of the particle size distribution function as the cloud microphysics scheme, including indirect effects, and then passed to the mstrnX scheme. In contrast, fixed effective radii of 8 µm for liquid hydrometeors and 40 µm for ice hydrometeors were assumed in NICAM.12. The use of consistent assumptions of coupling between cloud microphysics and radiative transfer can reduce the model bias in the radiation budget (Seiki et al., 2015a) and has nonnegligible impacts on climate projection (Chen et al., 2016). In addition, the coupling provides model developers with a better understanding of the origins of model biases (Hashino et al., 2016).

The mstrnX scheme requires a database of single-scattering properties of hydrometeors (RADPARA), including parameters such as the volume extinction coefficient, absorption coefficient, asymmetry factor, and second moment of phase function (Nakajima et al., 2000). In NICAM16-S we used the RADPARA database revised by Seiki et al. (2014). The RADPARA database of liquid hydrometeors was precalculated according to the Mie theory. The nonspherical RADPARA database developed by Fu (1996) and Fu et al. (1998) was applied to solid hydrometeors. The RADPARA database was then compiled as a lookup table of the effective radii from 1 µm to 1 mm to cover the size range of most of the hydrometeors in global simulations (Seiki et al., 2014). The effects of precipitating hydrometeors on the radiation budget are detectable, specifically over the intertropical convergence zone (ITCZ) and storm-track region (e.g., Waliser et al., 2011; Li et al., 2014, 2016; Chen et al., 2018; Michibata et al., 2019). The revised RADPARA database was evaluated in depth by comparing it with balloon-borne sonde observations in a midlatitude cirrus case (Seiki et al., 2014), and its effectiveness for global simulations was evaluated in several studies (Seiki et al., 2015a, b; Satoh et al., 2018).

Because of nonsphericity, the effective radius of ice particles has a controversial definition, whereas the effective radius is well defined in the case of spherical particles. According to Fu (1996), the effective radius, r_e, of solid hydrometeors is defined as follows:

$$\begin{array}{}\text{(1)}& r_{e,j}={\displaystyle \frac{\mathrm{3}}{\mathrm{4}{\mathit{\rho }}_{\mathrm{ice}}}}{\displaystyle \frac{\mathit{\rho }{q}_j}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{A}_j\left(D_j\right)N_j\left(D_j\right)\mathrm{d}{D}_j}}\left(j=i,s,g\right),\end{array}$$

where $j=i,s,g$ are cloud ice, snow, and graupel, respectively; ρ is the air density; ρ_ice=916.7 kg m⁻³; q_j is specific content; A_j is the projected area of a particle to flow; and D_j is diameter. The integral in Eq. (1) is analytically calculated using the assumed particle size distribution functions N_j(D_j) and a sponge-like spherical shape for cloud ice and graupel. In contrast, snow has two-dimensional fractal shapes; hence, the numerator–denominator ratio becomes almost constant. Thus, the effective radius of snow is assumed to be constant ($r_{\mathrm{e},\mathrm{s}}=\mathrm{125}$ µm with $A_{\mathrm{s}}=\mathrm{0.45}{D}_{\mathrm{s}}^{\mathrm{2.0}})$, and it is derived by approximating the A–D relationship of aggregates compiled by Mitchell (1996).

The impacts of the coupling procedure and the nonspherical scattering were examined by comparing the REF run with the NONSI run using the fixed effective radii and the spherical RADPARA database. The impacts of nonspherical scattering alone could be seen from the comparison between the REF and NONS runs using the spherical RADPARA database. Figure 8 shows the zonal mean values of the OLR and the outgoing shortwave radiation (OSR) at the TOA from the sensitivity experiments. Given the substantial increase in cloud ice and snow from the revised NSW6 scheme (cf. Sect. 3.2), ice optical thickness increases proportionately to the increases in the ice water path with fixed effective radii. As a result, both the longwave and shortwave radiation budgets are strongly biased in the NONSI run. A major portion of the biases in OLR is drastically offset in the NONS run by assuming larger effective radii in the coupling procedure. Thus, the coupling procedure automatically prevents artificial biases originating from inconsistent parameter settings between the cloud microphysics and radiative transfer with the model update. The use of the nonspherical RADPARA database slightly increases OSR over the tropics to the midlatitude because the assumed asymmetry factor for nonspherical particles is smaller than that for spherical particles (cf. Seiki et al., 2014).

Figure 8 Annual mean of the zonal mean OLR (a) and OSR (b) at the TOA (W m⁻²) for the CERES product (black) and NICAM16-9S. Green, blue, and red curves show the NONSI, NONS, and REF runs, respectively.

Finally, NICAM16-S still shows strong negative biases in OLR over the tropical to subtropical regions and OSR over the subtropical high-pressure belt at the TOA compared with the Clouds and Earth's Radiant Energy Systems (CERES) product (Fig. 8, black versus red curves), and these biases are qualitatively similar to those simulated in NICAM.12 (Kodama et al., 2015). The former OLR bias can be solved by increasing vertical resolution to 400 m near the tropopause with 74 vertical layers (Seiki et al., 2015b). The latter OSR bias mainly stems from the underestimation of low-level clouds, and the current updates in the cloud microphysics scheme do not improve the results for warm clouds. The coupling procedure strengthens the negative biases in OSR at TOA in the midlatitudes (Fig. 8, green versus blue curves) because, in NICAM.12, high clouds associated with extratropical cyclones are artificially brightened and, therefore, conceal the biases due to low-level clouds (cf. Fig. 4 in Kodama et al., 2012). Unlike NICAM.12, a strong negative bias of OSR at the TOA is also prominent over the Arctic region, and this seems to relate to the update (reduction) of the surface albedo introduced in Sect. 3.6.

3.4 Aerosols in the cloud microphysics and radiation schemes

In NICAM.12, the direct radiative effect of aerosols is not considered in the radiation scheme, and the number concentration of CCN is set to a constant value of 50 cm⁻³, a typical value over the ocean, in the cloud microphysics scheme. In NICAM16-S, both the direct and indirect effects of aerosols are considered by prescribing a distribution of aerosol mass concentration in the radiation scheme and CCN in the cloud microphysics scheme. The dataset of natural aerosols in the troposphere and stratosphere is described in Sect. 2.3. For anthropogenic aerosols in the troposphere, a simple plume model, MACv2-SP (Stevens et al., 2017; Fiedler et al., 2019), is used to diagnose the vertical profile of the aerosol optical depth, single-scattering albedo, asymmetry factor, and the factor of CCN increase arising from anthropogenic aerosols. This means that the magnitude of the anthropogenic increase in CCN depends on CCN from natural origins.

The sensitivity experiments with and without the update of the aerosol treatment (the REF and NOAER runs) are compared in Fig. 9. The aerosol update reduces the net downward shortwave radiation at the surface over most of the continents, particularly over Africa and South Asia (Fig. 9b and d), leading to a reduction of the excess insolation there (not shown). The reduced insolation at the surface is partly canceled out by a reduction in the upward longwave radiation over Africa (Fig. 9a and c) in association with a decrease in surface air temperature (not shown). Meanwhile, an enhancement of the surface net downward shortwave radiation is dominant over the ocean, particularly in the NICAM16-9S (14 km mesh) run (Fig. 9b and d), in association with a thinning of cloud optical depth and a decrease in cloud amount (not shown). This links to a decrease in the liquid water path with the aerosol update (Fig. 5f), which was also found in an online aerosol experiment with the 14 km mesh NICAM (Sato et al., 2018). As a result of these compensations, the global mean net surface radiation change arising from aerosol forcing is around −2.0 W m⁻² in NICAM16-7S and +0.4 W m⁻² in NICAM16-9S in this study. Such a reverse from a decrease to an increase between the resolutions might be related to the resolution dependency of the low and middle cloud amount in the REF run (left panels of Fig. 5h and i), and a detailed analysis is needed to properly understand the mechanism.

Figure 9 Impact of the aerosol treatment update (the REF run minus the NOAER run). Two-dimensional distribution, zonal, and global means of the net longwave (a, c) and shortwave (b, d) radiation at the surface are shown for the NICAM16-9S (a, b) and NICAM16-7S (c, d) runs. The data are gridded at 1^∘ in longitude and latitude, and the sign of the radiation is downward positive.

3.5 Land surface model

The land surface model MATSIRO (Takata et al., 2003) is used in NICAM. Recently, a wetland scheme was implemented in MATSIRO of NICAM16-S to represent the storage of snowmelt while considering the subgrid-scale terrain complexity (Nitta et al., 2017). The wetland scheme reduces the summertime warm and dry bias over much of Western Eurasia and North America through delayed snowmelt runoff in MIROC5 (Nitta et al., 2017). In addition, the effect of decreased surface albedo associated with the accumulation of water on land ice was implemented in NICAM16-S.

Figure 10 shows the impact of the land surface model update on soil moisture during boreal summers. The soil moisture is increased over most of the Eurasian and North American continents, as expected from Nitta et al. (2017), particularly in Siberia and around the Great Lakes. Though it is expected from Nitta et al. (2017) that increased soil moisture leads to an increase in precipitation and a decrease in surface air temperature, the actual impacts on precipitation and surface air temperature are still unclear (not shown). It is difficult to show robust reduction of the biases at this stage, and longer integration is needed to assess these impacts appropriately.

Figure 10 Impact of the land surface model update (the REF run minus the NOLND run) on the simulated soil moisture at the uppermost land surface model level. Simulations using NICAM16-7S were performed for 4 years, and June–July–August data for the last 3 years are averaged. The data are gridded at 2.5^∘ in longitude and latitude.

3.6 Surface albedo

Surface albedo values were updated based on observations. In NICAM.12, they were tuned to reduce the TOA radiation imbalance, and this caused a higher bias of surface albedo over the Arctic compared with that seen in satellite observations (Hashino et al., 2016). Table 6 shows the surface albedo values used in NICAM16-S and NICAM.12. The albedo values of the sea ice for the visible and near-infrared wavelengths and of fresh snow over land for the visible wavelengths were set to be smaller in NICAM16-S than those in NICAM.12. In addition to the updates in Table 6, an artificial elevation of the ocean surface albedo for the direct visible wavelength by a factor of 1.35 for the radiation scheme in NICAM.12 was discarded in NICAM16-S.

Table 6 Surface albedo in NICAM16-S and NICAM.12.

VIS, NIR, and IR stand for visible, near-infrared, and infrared bands, respectively.

The sensitivity experiments with and without the albedo update (the REF and NOALB runs, respectively) show that the use of the new surface albedo values tends to reduce the surface air temperature bias over land ice compared with the previous bias (Fig. 11a versus Fig. 11c and green versus red curves in Fig. 11d). Specifically, the cold bias in Greenland, the Himalayas, and the Antarctic is reduced. This is consistent with the reduced surface albedo for the visible wavelengths and the resulting decreased net upward shortwave radiation at the surface (Fig. 5k). Global mean net downward longwave radiation at the surface is increased (Fig. 5j). The increase is attributed primarily to a decrease in upward longwave radiation over the ocean (not shown), consistent with the increased surface albedo for the infrared wavelengths over the ocean. In terms of the TOA radiation budget, OSR worsens by a few watts per square meter (Fig. 5d), which arises from the polar regions (Fig. 12b; green versus red curves).

Figure 11 Bias in the simulated surface air temperature (K) by NICAM16-7S against JRA-55 reanalysis averaged for June 2004–May 2005. The NOALB (a), NOSIC (b), and REF (c) runs are shown. The zonal mean biases of the surface air temperature for the NOALB (green), NOSIC (blue), and REF (red) runs are shown in (d).

Figure 12 Same as Fig. 8 but for the NOALB (green), NOSIC (blue), and REF (red) runs by NICAM16-7S.

3.7 Treatment of oceans

A mixed-layer slab ocean model similar to that in McFarlane et al. (1992) was implemented in NICAM. The model predicts SST, ICE, snow over sea ice, and snow temperature by solving the heat balance between the ocean, sea ice, snow, and atmosphere. The depth of the slab ocean model is set to 15 m globally, considering the better performance of the simulated precipitation pattern (Kodama et al., 2015) and MJO (Grabowski, 2006). A simple nudging technique is used to force the predicted SST and ICE toward their reference states with a relaxation time of τ_SST and τ_ICE, respectively. Specifically, τ_SST=7 d and τ_ICE=0 (i.e., ICE was fixed to the boundary condition) were used in the slab ocean experiments of this study and in the previous climate simulation with NICAM.12 (Kodama et al., 2015). Both τ_SST and τ_ICE were set to zero in the fixed SST/ICE experiments, including the HighResMIP simulations.

In the slab ocean model implemented in NICAM.12 and NICAM16-S, SIC is diagnosed from ICE as

$$\begin{array}{}\text{(2)}& \mathrm{SIC}=\left\{\begin{array}{ll}\sqrt{\frac{\mathrm{ICE}}{\mathrm{SICCRT}}}& \text{for ICE}<\text{SICCRT},\\ \mathrm{1}& \text{for ICE}\ge \text{SICCRT},\end{array}\right.\end{array}$$

where SICCRT is a parameter in kg m⁻².

In many cases, including the HighResMIP protocol, only the SST and SIC data are provided to run the model, and, therefore, the ICE data prescribed for the model should be diagnosed from SIC data. In NICAM.12 and NICAM16-S, ICE is diagnosed simply as

$$\begin{array}{}\text{(3)}& \mathrm{ICE}=\mathrm{SICCRT}\times {\mathrm{SIC}}^{\mathrm{2}}.\end{array}$$

In the previous study using NICAM.12, we often set the value of SICCRT to 300 kg m⁻², considering Eq. (2). However, this situation leads to an underestimation of ICE over most of the sea ice areas and causes a warm bias over the Arctic (Kodama et al., 2015; Fig. 11b; blue curve in Fig. 11d). Based on an ocean model result (Hiroaki Tatebe, personal communication, 2018), we performed a series of preliminary annual-scale experiments using NICAM16-7S, with SICCRT values of 1600 and 3200, to improve the surface air temperature over the Arctic. As a result of this crude tuning, SICCRT is set to 1600 kg m⁻² in NICAM16-S. This led to a significant reduction in the warm bias (Fig. 11b, versus Fig. 11c; blue versus red curves in Fig. 11d) and an excess of OLR at TOA (blue versus red curves in Fig. 12a) over the Arctic.

Figure 13 Annual mean precipitation during June 2004–May 2005 for the GPCP product (a), the REFFIX runs (b, d, f), and the REFSLB runs (c, e, g). Results from NICAM16-7S (b, c), NICAM16-8S (d, e), and NICAM16-9S (f, g) are shown, and tropical mean precipitation is noted at the top-right corner of each panel.

In the HighResMIP protocol, SST and SIC in the model are fixed to the time-varying boundary conditions. Overall, the global mean impact of the slab ocean model is not very large (circles versus rectangles in Fig. 5). Compared with the fixed SST runs, global mean precipitation and OLR showed slight increases in the slab ocean runs, which are associated with a slightly warmer surface air temperature. In terms of the local climate, however, the fixed SST simulation is known to cause severe bias in the horizontal distribution of clouds and precipitation systems in the tropics (Kodama et al., 2015; Fig. 13), and thus the use of the slab ocean model with a 7 d relaxation time is often preferred. The introduction of the slab ocean model considerably affected the horizontal distribution of clouds and precipitation systems (Fig. 13). A double ITCZ bias was more prominent in the precipitation as well as the high cloud fraction and OLR fields (not shown) in the fixed SST runs compared with the slab ocean runs, particularly in the high-resolution run. Our investigation showed that NICAM16-9S with the slab ocean model best simulated the ITCZ peak precipitation and the precipitation pattern. Although the importance of the short-term SST variation driven by the atmosphere for the precipitation pattern is apparent from Fig. 13, the introduction of the slab ocean model alone does not resolve its bias. Further analysis is necessary to understand the physical mechanisms of the bias, which may be due to factors such as the convective timescale.

3.8 Orographic gravity wave drag

No gravity wave drag scheme is used in NICAM.12. In NICAM16-S, the conventional orographic gravity wave drag scheme (McFarlane, 1987) is used to better simulate the location and strength of the subtropical jet. The wave generation parameter α, which is proportional to the product of wave generation efficiency and representative horizontal wavenumber (Eq. 3.1b in McFarlane, 1987), was tuned first for NICAM16-9S to improve zonal mean zonal wind and then roughly halved as the horizontal mesh size was doubled. Specifically, α was set to $\mathrm{3.38}\times {\mathrm{10}}^{-\mathrm{5}}$ for NICAM16-7S, $\mathrm{7.12}\times {\mathrm{10}}^{-\mathrm{5}}$ for NICAM16-8S, and $\mathrm{1.46}\times {\mathrm{10}}^{-\mathrm{4}}$ for NICAM16-9S. Figures 14 and 15 show the zonal mean zonal wind in boreal summer and winter, as simulated with and without the gravity wave drag scheme (the REF and NOGWD runs, respectively). The zonal mean zonal wind simulated without the gravity wave drag scheme was biased poleward in both the NICAM16-9S and NICAM16-7S runs. The gravity wave drag scheme decelerated the zonal mean zonal wind at the poleward flank of the subtropical jet, especially during northern hemisphere winter, reducing the locational bias of the jet. The impact of the gravity wave drag scheme is larger in NICAM16-7S than in NICAM16-9S. The pattern of the zonal wind response to the gravity wave drag scheme is similar to that of previous studies (e.g., McFarlane, 1987; Iwasaki et al., 1989).

Figure 14 Zonal mean zonal wind (contours) and its bias from JRA-55 reanalysis (shadings) for June–August 2004 (a, b) and December 2004–February 2005 as simulated by NICAM16-9S without (a, c; NOGWD run) and with (b, d; the REF run) the gravity wave drag scheme.

Figure 15 Same as Fig. 14 but for the simulation by NICAM16-7S.

Although it is believed that even a mesh size of 14 km is insufficient to explicitly simulate the effects of the orographic gravity wave drag on the mean field (Nappo, 2012), introducing such a gravity wave drag scheme will not necessarily lead to an improvement of the simulated climate in the global nonhydrostatic model. The gravity wave drag scheme introduces the parameter α, which is uncertain and tuned to best simulate the climatology of the zonal wind for each resolution in general, although we did not tune α for each resolution. There is no solid guideline for determining α, including the dependency of the wave generation efficiency and representative horizontal wavenumber on the horizontal and vertical resolutions. Therefore, the use of a gravity wave drag scheme may hinder the pure resolution dependency of the mean field and suppress the advantages of the high-resolution model for simulating large-scale circulation in a seamless manner. Nevertheless, we decided to use the orographic gravity wave drag scheme for HighResMIP to reduce the locational bias of the subtropical jet to improve results for the tropical cyclone track region. It is important to recognize the merits and demerits involved in the use of a gravity wave drag scheme and reconsider its use depending on the main purpose of the simulation.

4 Horizontal and temporal resolution dependency

Understanding the dependency of horizontal resolution is a central interest of the HighResMIP. Figure 16 shows the global mean climate in NICAM16-7S (56 km mesh; blue circle), NICAM16-8S (28 km mesh; green circle), and NICAM16-9S (14 km mesh; red circle), along with its sensitivity to the time step of the models' dynamics (including the gravity wave drag scheme) and radiation scheme. Note that the dependency of the time step of the radiation scheme on the global mean climate is negligible.

Figure 16 Same as the left part of Fig. 5 but for the REFFIX, DDT2M, DDT1M, RDT20M, and RDT10M runs.

Global mean precipitation and TOA OLR decreased as the horizontal resolution increased (Fig. 16b and c), consistent with a previous NICAM study using 3.5 to 14 km meshes (Miyakawa and Miura, 2019). The results did not strongly depend on the time step of the model. As seen in Fig. 13, the precipitation pattern in the tropics depends strongly on the resolution: more-dominant double-ITCZ patterns and less-intense local precipitation are simulated as the horizontal resolution is increased. The intense precipitation occurs less frequently in the higher-resolution runs (Fig. 17), consistent with Noda et al. (2012) using older NICAM with 14 and 7 km meshes. The intense precipitation occurs more frequently in the model compared with the Global Precipitation Climatology Project (GPCP) product (Noda et al., 2012; Fig. 17), and it is consistent with the results of Maher et al. (2018), who compared precipitation in global climate models without a convection scheme with precipitation in the GPCP product. Na et al. (2020) showed that the frequency of intense precipitation in the GPCP product is lower than that in the TRMM product, and that the 14 km mesh NICAM without a convection scheme could realistically reproduce the intense precipitation observed by the TRMM satellite (Fig. 17).

Figure 17 Frequency of occurrence (%) of daily mean precipitation binned with an interval of 1 mm d⁻¹ during 1 June 2004–31 May 2005 as averaged over 15^∘ S–15^∘ N. The GPCP and TRMM products are shown with black and gray curves, respectively. The REFFIX runs with NICAM16-7S, NICAM16-8S, and NICAM16-9S are shown with blue, green, and red curves, respectively. The data were regridded at 1^∘ in longitude and latitude before the sampling.

The low cloud amount was substantially underestimated, especially in the NICAM16-7S and NICAM16-8S runs (Fig. 16i), leading to an underestimation of the OSR at the TOA (Fig. 16d). Consistently, the net downward shortwave radiation at the surface decreased (Fig. 16k) and the surface air temperature slightly decreased (Fig. 16a) as the horizontal resolution increased. This dependency of the low cloud amount and its related variables, in terms of global mean, on horizontal resolution could be reproduced by changing the time step of the dynamics in the model (the DDT2M and DDT1M runs). Specifically, the time loop in the model is based on the dynamics (Sect. 2.1), and the number of subcycles of the microphysics, turbulence, and surface processes were reduced in the DDT2M and DDT1M runs to keep their time steps unchanged (Table 2). This leads to an increase in frequency of the coupling between the dynamics and these physics schemes. Additional monthly scale sensitivity experiments suggest that the differences between the REFFIX runs and the DDT2M and DDT1M runs are mostly attributed to the altered frequency of physics–dynamics coupling, not the stability issue of the dynamical core (not shown). The low cloud amount was rather greater in the 56 km mesh run than that in the 14 km mesh run under the fixed time step of the dynamics (red circle in the REFFIX run versus blue circle in the DDT1M run in Fig. 16i). Also, the simulated OSR at the TOA is greater and closer to that of the CERES product in the 56 km mesh run compared with the 14 km mesh run with the same temporal resolution, though the better performance in the simulated global mean OSR at the TOA in the 56 km mesh run is a result of a strong compensation between a negative bias off the subtropical west coasts of continents and the southern hemisphere storm-track region and a positive bias in the rest of the lower latitudes (not shown). Such a result of dependency of OSR at the TOA on horizontal resolution with a fixed temporal resolution is similar to that found by Goto et al. (2020), who executed 14 and 56 km mesh online-aerosol NICAM with the same time steps of 60 s for the dynamics, turbulence, and surface schemes and 10 s for the cloud microphysics scheme. Beside the resolution, HighResMIP requires that the model parameters and configuration be the same between the standard and high-resolution runs. However, adapting the time step to the horizontal resolution, as we have done, is a common approach mainly for numerical stability reasons. Therefore, we highlighted here the impact of both the horizontal and temporal resolutions on the global mean climate in the NICAM simulations.

5 Computational aspects

5.1 Simulations

Table 7 describes the computational setting and the simulated year per wall-clock day (SYPD) of the NICAM16-S simulations run on the Earth Simulator 3. The Earth Simulator 3 is an NEC SX-ACE system with 5120 nodes in total for computation, and each computation node has four cores. We often used 10, 40, and 160 computation nodes to run NICAM16-7S, NICAM16-8S, and NICAM-9S, respectively, to obtain a balance between computational efficiency and wall-clock time. An exception was NICAM16-8S for the HighResMIP simulation, for which 160 computation nodes were used to perform a 101-year simulation within a realistic clock time. A file-staging option was used in the NICAM16-8S and NICAM16-9S simulations, whereas the file was directly read and written from a global file system in the NICAM16-7S simulation to reduce queuing time. The actual SYPD was a few times smaller than that shown in Table 7 for NICAM16-8S and NICAM16-9S.

Table 7 Computational aspects of the simulations on the Earth Simulator 3. They are sampled from 6-month simulations (1 July 2004–31 December 2004).

Figure 18 breaks down the total elapsed time for each component of the model. Unlike a previous evaluation using the K computer (Yashiro et al., 2016), the measurement includes the initial setup and input/output processes. Physics contributes a major part to the total elapsed time. Among the several physics components, the radiation scheme primarily contributes to the total elapsed time, followed by the cloud microphysics scheme, consistent with Yashiro et al. (2016) on the K computer. As the resolution increased, the percentage of time consumed by dynamics increased and the time consumed by the cloud microphysics and turbulence processes decreased because of their invariant time step among the models with different resolutions. An increase in the percentage of the land surface scheme in the simulations with higher resolution and greater number of computational nodes seems to be caused by the node imbalance associated with land–ocean distribution. Because the dynamics involves communication at every time step, some parts of the elapsed time counted as the dynamics can actually be a waiting time due to the node imbalance occurring in other processes.

Figure 18 Percentage of elapsed time for each component of NICAM16-S on the Earth Simulator 3. Times were sampled from 6-month simulations (1 July 2004–31 December 2004). The computational time per 1 d integration is shown on the left.

5.2 Postprocessing

The data are output on the model's native icosahedral grid on the height above sea level (ASL) or on the standard pressure. The vertical interpolation from the terrain-following height to the ASL height or standard pressure is performed online using second-order Lagrange interpolation during the simulation.

Postprocessing is performed in the following order: ico2ll, roughen, and z2pre.

1. ico2ll

   All the native icosahedral grid data are converted to high-resolution latitude–longitude grid data by area–weight averaging. The interval of latitude and longitude is determined so that the longitudinal interval is close to the average interval of the icosahedral grid (Satoh et al., 2014) on the Equator. Specifically, the interval of longitude and latitude is 0.56^∘ for NICAM16-7S, 0.28^∘ for NICAM16-8S, and 0.14^∘ for NICAM16-9S.

2. roughen

   All the high-resolution latitude–longitude grid data are coarsened to low-resolution latitude–longitude grid data by area–weight averaging. It is often necessary to reduce the amount of data by such coarsening, although HighResMIP does not ask us to coarsen the data. We prepared 1.0, 1.25, and 2.5^∘ data for analysis.

3. z2pre

   Several three-dimensional variables are converted from the ASL height to the standard pressure at this point by linear interpolation, if it is necessary.

After postprocessing steps (1)–(3), monthly mean data are created.

The pressure velocity is obtained from the vertical velocity and temperature after ico2ll under the assumption of hydrostatic balance. The geopotential height is calculated from the linear interpolation of vertical levels using logarithms of pressure, with an assumption of constant gravity acceleration regardless of height, which is a treatment consistent with the model configuration.

6 Summary

This paper describes the experimental design, the latest stable version of NICAM prepared for the HighResMIP (NICAM16-S), and impacts of NICAM updates on the simulated climatology using NICAM16-S at different resolutions. The major updates and their impacts on simulation results are summarized as follows:

• Update of the cloud microphysics scheme: snow and cloud ice were increased, leading to decreased high cloud amount and increased OLR at the TOA.

• Implementation of the coupling between cloud microphysics and radiation schemes: the negative OLR bias was reduced in association with a larger cloud ice effective radius.

• Update of treatment of natural and anthropogenic aerosols: the local surface radiation budget was improved, especially over Africa and South Asia.

• Update of land surface model: overall, the soil moisture over most of the Eurasian and the North American continents was increased.

• Update of the surface albedo values: the cold bias in Greenland, the Himalayas, and the Antarctic was decreased.

• Modification of the ICE diagnostics: the warm bias over the Arctic region was decreased.

• Introduction of the orographic gravity wave drag scheme: the location and strength of the zonal mean jet were improved.

Comprehensive evaluations and future projections using full HighResMIP data by NICAM16-S will be presented in a forthcoming paper.