Description and basic evaluation of simulated mean state, internal variability, and climate sensitivity in MIROC6

The sixth version of the Model for Interdisciplinary Research on Climate (MIROC), called MIROC6, was cooperatively developed by a Japanese modeling community. In the present paper, simulated mean climate, internal climate variability, and climate sensitivity in MIROC6 are evaluated and briefly summarized in comparison with the previous version of our climate model (MIROC5) and observations. The results show that the overall reproducibility of mean climate and internal climate variability in MIROC6 is better than that in MIROC5. The tropical climate systems (e.g., summertime precipitation in the western Pacific and the eastward-propagating Madden–Julian oscillation) and the midlatitude atmospheric circulation (e.g., the westerlies, the polar night jet, and troposphere–stratosphere interactions) are significantly improved in MIROC6. These improvements can be attributed to the newly implemented parameterization for shallow convective processes and to the inclusion of the stratosphere. While there are significant differences in climates and variabilities between the two models, the effective climate sensitivity of 2.6 K remains the same because the differences in radiative forcing and climate feedback tend to offset each other. With an aim towards contributing to the sixth phase of the Coupled Model Intercomparison Project, designated simulations tackling a wide range of climate science issues, as well as seasonal to decadal climate predictions and future climate projections, are currently ongoing using MIROC6.

1 Introduction

As global warming due to increasing emissions of anthropogenic greenhouse gases progresses, global and regional patterns of atmospheric circulation and precipitation as well as temperature are projected to be drastically changed by the end of the twenty-first century (e.g., Neelin et al., 2006; Zhang et al., 2007; Bengtsson et al., 2009; Andrews et al., 2010; Scaife et al., 2012); the occurrence frequency of extreme weather events such as heat waves and droughts will be increased, and extratropical cyclones will be stronger than in the present (e.g., Mizuta, 2012; Sillmann et al., 2013; Zappa et al., 2013). Corresponding to the atmospheric changes under global warming, the sea levels will rise due to the thermal expansion of seawater and ice sheet melting in the polar continental regions (e.g., Church and White, 2011; Bamber and Aspinall, 2013). Additionally, ocean acidification due to the absorption of atmospheric carbon dioxide (CO₂) and changes in carbon–nitrogen cycles are expected to lead to the loss of Earth biodiversity (e.g., Riebesell et al., 2009; Rockström et al., 2009; Taucher and Oschlies, 2011; Watanabe and Kawamiya, 2017). Societal demands for information on global and regional climate changes have increased significantly worldwide in order to meet information requirements for political decision-making related to mitigation and adaptation to global warming.

The Intergovernmental Panel on Climate Change (IPCC) has continuously published assessment reports (ARs) in which a comprehensive view of past, present, and future climate changes on various timescales, including centennial global warming, is synthesized. Together with observations, climate models have been contributing to the IPCC ARs through a broad range of numerical simulations, especially future climate projections after the twenty-first century. However, there are many uncertainties in future climate projections, and the range of uncertainties has not been narrowed by an update of the IPCC reports. The uncertainties are arising from imperfections of climate models in representing microscale to global-scale physical and dynamical processes in subsystems of the Earth's climate and their interactions. To reduce the uncertainties and errors in climate projections and predictions, utilizing observations, extracting the essences of physical processes in the real climate, and investigating the response of the climate system to various external forcings based on a set of climate model simulations are necessary. In particular, a state-of-the-art climate model that can represent various processes in the Earth's climate system is a powerful tool for a deeper understanding of the Earth's climate system.

One Japanese climate model, which is called MIROC (Model for Interdisciplinary Research on Climate), has been cooperatively developed at the Center for Climate System Research (CCSR; the precursor of a part of the Atmosphere and Ocean Research Institute), the University of Tokyo, the Japan Agency for Marine-Earth Science and Technology (JAMSTEC), and the National Institute for Environmental Studies (NIES). Utilizing MIROC, our Japanese climate modeling group has been tackling a wide range of climate science issues as well as seasonal to decadal climate predictions and future climate projections. At the same time, by providing simulation data, we have been participating in the third and fifth phases of the Coupled Model Intercomparison Project (CMIP3 and CMIP5; Meehl et al., 2007; Taylor et al., 2012) that has been contributing to the IPCC ARs by synthesizing multi-model ensemble datasets.

In the years up to the IPCC fifth assessment report (IPCC AR5; IPCC, 2013), we have developed four versions of MIROC, three of which (MIROC3m, MIROC3h, and MIROC4h) have almost the same dynamical and physical packages but different resolutions. MIROC3m (K-1 model developers, 2004) is composed of T42L20 atmosphere and 1.4^∘L43 ocean. Resolutions of MIROC3h (K-1 model developers, 2004) are higher than MIROC3m and are T106L56 for the atmosphere and eddy-permitting for the ocean ($\mathrm{1}/\mathrm{4}{}^{\circ }\times \mathrm{1}/\mathrm{6}{}^{\circ }$). Only the horizontal resolution of the atmosphere of MIROC3h is changed to T213 in MIROC4h (Sakamoto et al., 2012). MIROC5 is composed of T85L40 atmosphere and 1.4^∘L50 ocean but with considerably updated physical and dynamical packages (Watanabe et al., 2010). These models have been used to study various scientific issues such as the detection of natural influences on climate changes (e.g., Nozawa et al., 2005; Mori et al, 2014; Watanabe et al., 2014), uncertainty quantification of climate sensitivity (e.g., Shiogama et al., 2012; Kamae et al., 2016), future projections of regional sea level rises (e.g., Suzuki et al., 2005; Suzuki and Ishii, 2011), and mechanism studies on tropical decadal variability (e.g., Tatebe et al., 2013; Mochizuki et al., 2016).

During the last decade, our efforts have been preferentially devoted to providing science-oriented risk information on climate changes that is beneficial to international, domestic, and municipal communities. For example, so-called event attribution (EA) studies with large-ensemble simulations initiated from slightly different conditions have been conducted in order to statistically evaluate the influences of global warming on the occurrence frequencies of observed individual extremes (e.g., Imada et al., 2013; Watanabe et al., 2013; Shiogama et al., 2014). Seasonal to decadal climate predictions are also of significant concern. By initializing prognostic variables in our climate models using observation-based data (Tatebe et al., 2012), significant prediction skills in several specific phenomena, such as the El Niño–Southern Oscillation (ENSO) and the Arctic sea ice extent on seasonal timescales, the Pacific Decadal Oscillation (PDO; Mantua et al., 1997), the Atlantic Multi-decadal Oscillation (AMO; Schlesinger and Ramankutty, 1994), and the tropical trans-basin interactions between the Pacific and the Atlantic on decadal timescales, are detected (e.g., Mochizuki et al., 2010; Chikamoto et al., 2015; Imada et al., 2015; Ono et al., 2018).

However, while the applicability of MIROC has been extended to a wide range of climate science issues, almost all of the abovementioned approaches were based on our medium-resolution versions of MIROC (MIROC3m and MIROC5), and it is well known that higher-resolution models are capable of better representing the model mean climate and internal climate variability, such as regional extremes, orographic winds, and oceanic western boundary currents and eddies, than lower-resolution models (e.g., Shaffrey et al., 2009; Roberts et al., 2009; Sakamoto et al., 2012). Nevertheless, persistent biases remain associated with, for example, cloud–aerosol–radiative feedback and turbulent vertical mixing of the air in the planetary boundary layer (e.g., Bony and Dufresne, 2005; Bodas-Salcedo et al., 2012; Williams et al., 2013), which are tightly linked with dominant uncertainties in climate projections. Therefore, improvement of physical parameterizations for sub-grid-scale processes is essential for better representing observed climatic mean states and internal climate variability. In addition to physical parameterizations, enhanced vertical resolution in both atmosphere and ocean components, along with a highly accurate tracer advection scheme, has been suggested to have impacts on the reproducibility of model climate and internal climate variations (e.g., Tatebe and Hasumi, 2010; Ineson and Scaife, 2009; Scaife et al., 2012).

Recently, we have developed the sixth version of MIROC, called MIROC6. This newly developed climate model has updated physical parameterizations in all sub-modules. In order to suppress an increase in computational cost, the horizontal resolutions of MIROC6 are not significantly higher than those of MIROC5. The reason is that a larger number of ensemble members is required to realize significant seasonal predictions of, for example, the wintertime Eurasian climate (Murphy, 1990; Scaife et al., 2014). Indeed, climate predictions by the older versions of MIROC having at most 10 ensemble members are skillful only in the tropical climate and the midlatitude ocean, not in the midlatitude atmosphere. Large-ensemble predictions are also required in decadal-scale predictions in order to evaluate the human influences on near-term climate changes. The model top in MIROC6 is placed at the 0.004 hPa pressure level, which is higher than that of MIROC5 (3 hPa), and the stratospheric vertical resolution has been enhanced in comparison to MIROC5 in order to represent the stratospheric circulation. Overall, the reproducibility of the mean climate and internal variability of MIROC6 is better than that of MIROC5, but the model's computational cost is about 3.6 times as large as that of MIROC5. Considering that the computational costs of large-ensemble predictions based on climate models with horizontal resolutions of, for example, 50 km atmosphere and eddy-resolving ocean are still huge on recent computer systems, the use of relatively low-resolution models such as MIROC6 with further elaborated parameterizations can still be actively useful in science-oriented climate studies and climate predictions produced for societal needs.

The rest of the present paper is organized as follows. We describe the model configuration, tuning, and spin-up procedures in Sect. 2, while simulated mean state, internal variability, and climate sensitivity are evaluated in Sect. 3. Simulation performance of MIROC6 and remaining issues are briefly summarized and discussed in Sect. 4. Currently, MIROC6 is being used for various simulations designed by the sixth phase of CMIP (CMIP6; Eyring et al., 2016), which aims to strengthen the scientific basis of IPCC AR6. Large-ensemble simulations and climate predictions using MIROC6 are also ongoing for science-oriented studies in our modeling group and for societal benefits. In addition, the latest Earth system model version of MIROC with the global carbon cycle, whose physical core will be MIROC6, has been developed for CMIP6 towards a further wide range of issues regarding climate and societal applications (Hajima et al., 2019).

2 Model configurations and spin-up procedures

MIROC6 is composed of three sub-models: atmosphere, land, and sea ice–ocean. The atmospheric model is based on the CCSR-NIES atmospheric general circulation model (AGCM; Numaguti et al., 1997). The land surface model is based on Minimal Advanced Treatments of Surface Interaction and Runoff (MATSIRO; Takata et al. 2003), which includes a river-routing model from Oki and Sud (1998) based on a kinematic wave flow equation (Ngo-Duc et al., 2007) and a lake module in which one-dimensional thermal diffusion and mass conservation are considered. The sea ice–ocean model is based on the CCSR Ocean Component model (COCO; Hasumi, 2006). A coupler system calculates heat and freshwater fluxes between the sub-models in order to ensure that all fluxes are conserved within machine precision and then exchanges the fluxes among the sub-models (Suzuki et al., 2009). No flux adjustments are used in MIROC6. In the remaining part of this section, we will provide details on MIROC6 configurations, focusing on updates from MIROC5. Readers may also refer to Table A1 in Appendix A where the updates are briefly summarized.

2.1 Atmospheric component

MIROC6 employs a spectral dynamical core in its AGCM component as in MIROC5. The horizontal resolution is a T85 spectral truncation that is an approximately 1.4^∘ grid interval for both latitude and longitude. The vertical grid coordinate is a hybrid σ–p coordinate (Arakawa and Konor, 1996). The model top is placed at 0.004 hPa, and there are 81 vertical levels (Fig. 1a). The vertical grid arrangement in MIROC6 is considerably enhanced in comparison to that in MIROC5 (40 levels; 3 hPa) so that the stratospheric circulation can be represented. A sponge layer that damps wave motions is set at the model-top level by increasing Rayleigh friction to prevent extra wave reflection near the model top. The atmospheric component of MIROC6 has standard physical parameterizations for cumulus convection, radiation transfer, cloud microphysics, turbulence, and gravity wave drag. It also has an aerosol module. These are basically the same as those used in MIROC5, but several updates have been made, as will be detailed below. The parameterizations for cloud microphysics and planetary boundary layer processes in MIROC6 are the same as in MIROC5. The standard time step for MIROC6 is 6 min, which is shorter than that of MIROC5 (12 min) because stratospheric winds whose speed sometimes exceeds 150 m s⁻¹ must be resolved in time integration. The time step for radiative transfer models is set separately and is 3 h in both MIROC6 and MIROC5.

Figure 1 Vertical half-levels for the atmospheric (a) and the oceanic (b) components of MIROC6 and MIROC5.

A cumulus parameterization proposed by Chikira and Sugiyama (2010), which uses an entrainment formulation from Gregory (2001), is adopted in MIROC6 as in MIROC5. This parameterization deals with multiple cloud types including shallow cumulus and deep convective clouds. MIROC5, however, tends to overestimate low-level cloud amounts over the low-latitude oceans and has a dry bias in the free troposphere. These biases appear to be the result of insufficient vertical mixing of the humid air in the planetary boundary layer and the dry air in the free troposphere. To alleviate these biases, an additional parameterization for shallow cumulus convection based on Park and Bretherton (2009) is implemented in MIROC6. Shallow convection associated with atmospheric instability is calculated by the Chikira and Sugiyama (2010) scheme, and that associated with turbulence in the planetary boundary layer is represented by the Park and Bretherton (2009) scheme. The shallow convective parameterization is a mass flux scheme based on a buoyancy-sorting, entrainment–detrainment, single-plume model that calculates the vertical transport of liquid water, potential temperature, total water mixing ratio, and horizontal winds in the lower troposphere. The cloud-base mass flux is controlled by turbulent kinetic energy within the sub-cloud layer and convective inhibition. The cloud-base height for shallow cumulus is set between the lifting condensation level and the boundary layer top, which is diagnosed based on the vertical profile of relative humidity. When implementing the parameterization in MIROC6, the following conditions for triggering the shallow convection are specified: (1) the estimated inversion strength (Wood and Bretherton, 2006) is smaller than a tuning parameter, and (2) the convection depth diagnosed by a separate cumulus convection scheme (Chikira and Sugiyama, 2010) is smaller than a tuning parameter.

The Spectral Radiation Transport Model for Aerosol Species (SPRINTARS; Takemura et al., 2000, 2005, 2009) is used as an aerosol module for MIROC6 to predict the mass mixing ratios of the main tropospheric aerosols, which are black carbon, organic matter, sulfate, soil dust, sea salt, and the precursor gases of sulfate (sulfur dioxide, SO₂, and dimethylsulfide). By coupling the radiation and cloud–precipitation schemes in MIROC, SPRINTARS calculates not only the aerosol transport processes of emission, advection, diffusion, sulfur chemistry, wet deposition, dry deposition, and gravitational settling, but also the aerosol–radiation and aerosol–cloud interactions. There are two primary updates in SPRINTARS of MIROC6 that were not included in MIROC5. One is the treatment of precursor gases of organic matter as prognostic variables. In the previous version, the conversion rates from the precursor gases (e.g., terpene and isoprene) to organic matter are prescribed (Takemura et al., 2000), while an explicit simplified scheme for secondary organic matter was introduced from a global chemical climate model (Sudo et al., 2002). The other is a treatment of oceanic primary and secondary organic matter. Emissions of primary organic matter are calculated with wind at a 10 m height, the particle diameter of sea salt aerosols, and chlorophyll a concentration at the ocean surface (Gantt et al., 2011). The oceanic isoprene and monoterpene, which are precursor gases of organic matter, are emitted depending on the photosynthetically active radiation, diffuse attenuation coefficient at 490 nm, and the ocean surface chlorophyll a concentration (Gantt et al., 2009).

The radiative transfer in MIROC6 is calculated by an updated version of the k-distribution scheme used in MIROC5 (Sekiguchi and Nakajima, 2008). The single-scattering parameters have been calculated and tabulated in advance, and liquid, ice, and five aerosol species can be treated in this updated version. Given the significant effect of crystal habit on a particle's optical characteristics (Baran, 2012), the assumption of ice particle habit has been updated from our previous simple assumption of a sphere used in MIROC5 to a hexagonal solid column (Yang et al., 2013) in MIROC6. The upper limits of the mode radius of cloud particles have been extended from 32 µm to 0.2 mm for liquids and from 80 µm to 0.5 mm for ice. Therefore, the scheme can now handle the large-sized water particles (e.g., drizzle and rain) that have been shown to have significant radiative impacts (Waliser et al., 2011).

Following Hines (1997) and Watanabe et al. (2011), a non-orographic gravity wave parameterization is newly implemented into MIROC6 in order to represent realistic large-scale circulation and thermal structures in the stratosphere and mesosphere. Following Watanabe (2008), a present-day climatological source of non-orographic gravity waves, which is estimated using the results of a gravity-wave-resolving version of the MIROC AGCM (Watanabe et al., 2008), is launched at the 70 hPa level in the extratropics, while an isotropic source of non-orographic gravity waves is launched at the 650 hPa level in the tropics. Together with this parameterization, an orographic gravity wave parameterization from McFarlane (1987) is also adopted as in MIROC5. In both the orographic and non-orographic gravity wave parameterizations, wave source parameters at launch levels are tuned so that the realistic seasonal progress of the middle atmosphere circulation, frequency of sudden stratospheric warmings, and period and amplitude of the equatorial quasi-biennial oscillations (QBOs) can be represented.

2.2 Land surface component

The land surface model is also basically the same as in MIROC5. Energy and water exchanges between land and atmosphere are calculated, considering the physical and physiological effects of vegetation with a single-layer canopy, as well as the thermal and hydrological effects of snow and soil, respectively, with three-layer snow and six-layer soil down to a 14 m depth. Sub-grid fractions of land use and snow cover have also been considered. The time step for the land surface model integration is 1 h in MIROC6, which is the same as in MIROC5. In addition to the standard package in MIROC5, a few other physical parameterizations are implemented as described below.

A physically based parameterization of sub-grid snow distribution (SSNOWD; Liston, 2004; Nitta et al., 2014) replaces the simple functional approach of snow water equivalent in calculating sub-grid snow fractions in MIROC5 in order to improve the seasonal cycle of snow cover. In SSNOWD, the snow cover fraction is formulated for accumulation and ablation seasons separately. For the ablation season, the snow cover fraction decreases based on the sub-grid distribution of the snow water equivalent. A lognormal distribution function is assumed and the coefficient of variation category is diagnosed from the standard deviation of the sub-grid topography, coldness index, and vegetation type that is a proxy for surface winds. While the cold degree month was adopted for coldness in the original SSNOWD, we decided instead to introduce the annually averaged temperature over the latest 30 years using the time relaxation method of Krinner et al. (2005), in which the timescale parameter is set to 16 years. The temperature threshold for a category diagnosis is set to 0 and 10 ^∘C. In addition, a scheme representing a snow-fed wetland that takes into consideration sub-grid terrain complexity (Nitta et al., 2017) is incorporated. The river-routing model and lake module are the same as those used in MIROC5, but the river network map is updated to keep the consistency with the new land–sea mask (Yamazaki et al., 2009).

2.3 Ocean and sea ice component

The ocean component of MIROC6 is basically the same as that used in MIROC5, but several updates are implemented as described below. The warped bipolar horizontal coordinate system in MIROC5 has been replaced by the tripolar coordinate system proposed by Murray (1996). Two singular points in the bipolar region to the north of about 63^∘ N are placed at (63^∘ N, 60^∘ E) in Canada and (63^∘ N, 120^∘ W) in Siberia (Fig. 2). In the spherical coordinate portion to the south of 63^∘ N, the longitudinal grid spacing is 1^∘ and the meridional grid spacing varies from about 0.5^∘ near the Equator to 1^∘ in the midlatitudes. In the central Arctic Ocean where the bipole coordinate system is applied, the grid spacings are about 60 km zonal and 33 km meridional, respectively. By introducing the horizontal tripolar coordinate system, it is expected that theoretical westward propagation of the oceanic baroclinic Rossby can be represented with fewer numerical dispersions because of agreement of the coordinate system and the geographical coordinate system. It is also expected that the horizontal resolutions in the Arctic Ocean where the Rossby radius of deformation is relatively small are higher than in the case in which the bipolar warped coordinate system in MIROC5 is adopted. There are 62 vertical levels in a hybrid σ–z coordinate system. The horizontal grid spacing in MIROC5 is nominally 1.4^∘, except for the equatorial region, and there are 49 vertical levels. The resolutions in MIROC6 are higher than in MIROC5. In particular, 31 (23) of the 62 (49) vertical layers in MIROC6 (MIROC5) are within the upper 500 m of depth (Fig. 1b). The increased number of vertical layers in MIROC6 has been adopted in order to better represent the equatorial thermocline and observed complex hydrography in the Arctic Ocean. An increase in computational costs of the ocean component due to higher resolutions in MIROC6 is suppressed by implementing a time-staggered scheme for the tracer and baroclinic momentum equations (Griffies et al., 2005). Owing to the time-staggered scheme, the time step for the ocean and sea ice components of MIROC6 is 20 min, which is longer than that in MIROC5 (15 min).

Figure 2 Horizontal grid coordinate system and model bathymetry of the ocean component of MIROC6.

The tracer advection scheme (Prather, 1986), the surface mixed layer parameterization (Noh and Kim, 1999), and the parameterization for eddy isopycnal diffusion (Gent et al., 1995) used in MIROC6 are the same as those used in MIROC5. Also as in MIROC5, the bottom boundary layer parameterization of Nakano and Suginohara (2002) is introduced south (north) of 54^∘ S (49^∘ N) to represent the downsloping flow of dense waters. The constant parameters used in the abovementioned parameterizations are determined in the same manner as that of MIROC5, except for the Arctic region. An empirical profile of background vertical diffusivity, which is proposed in Tsujino et al. (2000), is modified above the 50 m depth to the north of 65^∘ N. It is $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ in the uppermost 29 m and gradually increases to $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ at the 50 m depth. Additionally, the turbulent mixing process in the surface mixed layer is changed so that there is no surface wave breaking and no resultant near-surface mixing in regions covered by sea ice. The combination of the weak background vertical diffusivity and suppression of turbulent mixing under the sea ice contributes to better representations of the surface stratification in the Arctic Ocean with little impact on the rest of the global oceans (Komuro, 2014).

The sea ice component in MIROC6 is almost the same as in MIROC5. A brief description, along with some major parameters, is given here. Readers may refer to Komuro et al. (2012) and Komuro and Suzuki (2013) for further details. A sub-grid-scale sea ice thickness distribution is incorporated by following Bitz et al. (2001). There are five ice categories (plus one additional category for open water), and the lower bounds of the ice thickness for these categories are set to 0.3, 0.6, 1, 2.5, and 5 m. The momentum equation for sea ice dynamics is solved using elastic–viscous–plastic rheology (Hunke and Dukowicz, 1997). The strength of the ice per unit thickness and concentration is set at 2.0×10⁴ N m⁻², and the ice–ocean drag coefficient is set to 0.02. The surface albedo for bare ice surface is 0.85 (0.65) for the visible (infrared) radiation. The surface albedo in snow-covered areas is 0.95 (0.80) when the surface temperature is lower than −5 ^∘C for the visible (infrared) radiation, and it is 0.85 (0.65) when the temperature is 0 ^∘C. Note that the albedo changes linearly between −5 and 0 ^∘C. These parameter values listed here are the same as those listed in MIROC5.

2.4 Boundary conditions

A set of external forcing data recommended by the CMIP6 protocol is used. The historical solar irradiance spectra, greenhouse gas concentrations, anthropogenic aerosol emissions, and biomass burning emissions are given by Matthes et al. (2017), Meinshausen et al. (2017), Hoesly et al. (2018), and van Marle et al. (2017), respectively. The concentrations of greenhouse gases averaged globally and annually are given to MIROC6. Radiative forcing of stratospheric aerosols due to volcanic eruptions is computed by vertically integrating extinction coefficients for each radiation band, which are provided by Thomason et al. (2019), in the model layers above the tropopause. Three-dimensional atmospheric concentrations of historical ozone (O₃) are produced by the Chemistry–Climate Model Initiative (Hegglin et al., 2019; the data are available at http://blogs.reading.ac.uk/ccmi/forcing-databases-in-support-of-cmip6/, last access: 6 July 2016). Three-dimensional concentrations of the OH radical, hydrogen peroxide (H₂O₂), and nitrate (NO₃) are precalculated by a chemical atmospheric model from Sudo et al. (2002). As precursors of secondary organic aerosol, emission data on terpenes and isoprene provided by the Global Emissions Inventory Activity (Guenther et al., 1995) are normally used, although simulated emissions from the land ecosystem model of Ito and Inatmoni (2012) are also used alternatively.

For specifying the soil types and area fractions of natural vegetation and cropland on grids of the land surface component, the harmonized land use dataset (Hurtt et al., 2011), Center for Sustainability and the Global Environment global potential vegetation dataset (Ramankutty and Foley, 1999), and the dataset provided by the International Satellite Land Surface Climatology Project Initiative I (Sellers et al., 1996) are used. These datasets are also used in prescribing background reflectance at the land surface. Leaf area index data are prepared based on the moderate-resolution imaging spectroradiometer leaf area index products of Myneni et al. (2002).

The forcing dataset used for the preindustrial control simulation is basically composed of data for the year 1850, which are included in the abovementioned historical dataset. The stratospheric aerosols and solar irradiance in the preindustrial simulation are given as monthly climatology averaged in 1850–2014 and in 1850–1873, respectively. The total solar irradiance is about 1361 W m⁻², and the global mean concentrations of CO₂, methane (CH₄), and nitrous oxide (N₂O) are 284.32 ppm, 808.25 ppb, and 273.02 ppb, respectively.

2.5 Spin-up and tuning procedures

Firstly, the stand-alone ocean component of MIROC6, which includes the sea ice processes, is integrated from the initial motionless state with the observed temperature and salinity distribution of the Polar Science Center hydrographic climatology (Steele et al., 2001). Ocean model coastline geometry and bottom bathymetry are specified based on horizontal interpolation of the land and seafloor dataset of ETOPO5 (National Geophysical Data Center, 1993). The ocean component is spun up for 1000 years by the monthly climatological surface fluxes of Röske (2006). An acceleration method from Bryan (1984) is used in the spin-up stage in order to obtain a thermally and dynamically quasi-steady state. After the spin-up, additional integration for 200 years is performed without the acceleration method. By analyzing the last 50-year-long data from the stand-alone ocean component, the monthly climatology of typical variables (e.g., zonal mean temperature and salinity in several basins, volume transports across major straits and archipelagos, meridional overturning circulation, and sea ice distributions) is compared with observations. Once the configuration of the ocean component is frozen, the land–sea distribution and land–sea area ratios on the atmospheric and land surface model grids are determined according to the coastline geometry of the ocean component, after which the atmospheric and the land surface components are coupled with the ocean component. Surface topography in the atmospheric and land surface component is also made using the ETOPO5 dataset. Note that the horizontal grid arrangement of the land surface model is exactly the same as the atmospheric component. The coupling interval among the sub-models is 1 h. An initial condition of the ocean component in MIROC6 is given by the stand-alone ocean experiment, and those of the atmosphere and land are taken from an arbitrary year of the preindustrial control run of MIROC5.

After coupling the sub-models, climate model tuning is done under the preindustrial boundary conditions. Conventionally, the climate models of our modeling community are retuned in coupled modes after stand-alone sub-model tuning. This is because the reproducibility of climatic mean state and internal climate variations is not necessarily guaranteed in climate models with the same parameters determined in stand-alone sub-model tuning, which is particularly the case in the tropical climate. In our tuning procedures described below, many of the 10-year-long climate model runs are conducted with different parameter values. There are numerous parameters associated with physical parameterizations, whose upper–lower bounds are constrained by empirical or physical reasoning. The main parameters used in our tuning procedures are chosen by referring to a perturbed parameter ensemble set made by Shiogama et al. (2012) in which parameter sensitivity to cloud radiative processes is examined. The impact of parameter tuning on the present climate is also discussed by Ogura et al. (2017), focusing on top-of-atmosphere (TOA) radiation and clouds. Any objective and optimal methods for parameter tuning are not used in our modeling group, and the tuning procedures are like those in other climate modeling groups as summarized in Hourdin et al. (2017).

In the first model tuning step, climatology, seasonal progression, and internal climate variability in the tropical coupled system are tuned so that departures from observations or reanalysis datasets are reduced. Here, it should be noted that representation of the tropical system in MIROC6 is sensitive to the parameters for convection and planetary boundary layer processes. Specifically, parameters of reference height for cumulus precipitation, efficiency of the cumulus entrainment of the surrounding environment, and maximum cumulus updraft velocity at the cumulus base are used to tune the strength of the equatorial trade wind, the climatological position and intensity of the Intertropical Convergence Zone (ITCZ) and South Pacific Convergence Zone (SPCZ), and the interannual variability of the El Niño–Southern Oscillation (ENSO). In particular, the parameter for the cumulus entrainment is known as a controlling factor of ENSO in MIROC5 (Watanabe et al., 2011). Summertime precipitation in the western tropical Pacific that is characteristic of tropical intra-seasonal oscillations is tuned by using the parameter for shallow convection describing the partitioning of turbulent kinetic energy between horizontal and vertical motions at the sub-cloud layer inversion. Next, the wintertime midlatitude westerly jets and the stationary waves in the troposphere are tuned using the parameters of the orographic gravity wave drag and the hyper-diffusion of momentum. The parameters of the hyper-diffusion and the non-orographic gravity wave drag are also used when tuning stratospheric circulation of the polar vortex and QBO. Finally, the radiation budget at the TOA is tuned, primarily using the parameters for the auto-conversion process so that excess downward radiation can be minimized and maintained closer to 0.0 W m⁻². The surface albedos for bare sea ice and snow-covered sea ice are set to higher values than in observations (see Sect. 2.3) in order to avoid underestimating the summertime sea ice extent in the Arctic Ocean due to excess downward shortwave radiation in this region. In addition, parameter tuning for the total radiative forcing associated with aerosol–radiation and aerosol–cloud interactions is done. So that the total radiative forcing can be closer to the estimate of −0.9 W m⁻² (IPCC, 2013; negative value indicates cooling) with an uncertainty range of −1.9 to −0.1 W m⁻², parameters of cloud microphysics and the aerosol transport module, such as the timescale for cloud droplet nucleation, in-cloud properties of aerosol removal by precipitation, and the minimum threshold for the number concentration of cloud droplets, are perturbed. To determine a suitable parameter set, several pairs of a present-day run under the anthropogenic aerosol emissions at the year 2000 and a preindustrial run are conducted. A pair of present and preindustrial runs has exactly the same parameters, and differences of tropospheric radiation between two runs are considered anthropogenic radiative forcing. Note that MIROC6 in a coupled mode is used in this tuning procedure, and thus the sea surface temperature (SST) is not fixed. The estimated radiative forcing here is not strictly the same as the effective radiative forcing estimated in IPCC (2013). However, by the present tuning procedure, the global mean surface air temperature (SAT) change after the mid-nineteenth century is well reproduced in the historical runs by MIROC6 (details are discussed in Sect. 4). As mentioned above, the reproducibility of the global mean SAT is not a tuning goal but is a typical metric that reflects results of the parameter tunings for individual processes of convection, dynamics, and radiative forcing.

After fixing the model parameters, the climate model is spun up for 2000 years. During the first several hundred years, waters contained in the land surface are drained to the ocean via river runoff, which leads to a temporal weakening of the meridional overturning circulation in the ocean and a rising of the global mean sea level. After the global hydrological cycle reaches an equilibrium state, the strengths of the meridional overturning circulation recover and keep a quasi-steady state. The abovementioned processes take about 1000 years, after which an additional 1000-year-long integration is performed in order to obtain a thermally and dynamically quasi-steady ocean state.

Figure 3 shows the time series of the global mean quantities after the spin-up. The labeled year in Fig. 3 indicates the elapsed year after the spin-up duration of 2000 years. The linear trend of the global mean SAT is $\mathrm{9.5}\times {\mathrm{10}}^{-\mathrm{3}}$ K per century and is much smaller than the observed value of about 0.62 K per century in the twentieth century, indicating that there is no significant drift and the global mean SAT is in a quasi-steady state. While the global mean SST is in a quasi-steady state (linear trend of $\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{3}}$ K per century), the global mean ocean temperature shows a larger trend of $\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{3}}$ K per century in the first 500 years than that of $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$ K per century in the later period. In the later sections, the 200-year-long data between the 500th and 699th years are analyzed.

Figure 3 (a) Time series of the global mean SAT (solid) and the TOA radiation budget (dashed; upward positive). (b) Same as (a), but for the global mean SST (solid) and the ocean temperature through the full water column (dashed).

The trend of the global mean ocean temperature in the later period suggests slight but continuous warming of the deep ocean. The radiation budget at the TOA is 1.1 W m⁻² downward on average (linear trend of $\mathrm{9.5}\times {\mathrm{10}}^{-\mathrm{3}}$ K per century), and the net heat input at the sea surface is 0.32 W m⁻². The deep ocean warming is explained by the net heat input. Note that there is about 0.78 W m⁻² of inconsistency between the TOA radiation budget and the ocean heat uptake. This heat energy inconsistency is due to internal energy associated with precipitation, water vapor, and river runoff not being taken account in the atmospheric and land surface component in MIROC6, as well as the fact that these waters with no temperature information implicitly set their temperature to the SST when they flow or fall into the ocean. Perpetual melting of the prescribed Antarctic ice sheet with invariant ice thickness, which occurs due to the warm SAT bias in the Antarctic region (details will be discussed in Sect. 3.1.3), is also a cause of the heat energy inconsistency.

3 Results of preindustrial simulation

Representations of climatic mean field and internal climate variability in MIROC6 (Tatebe and Watanabe, 2018a) are evaluated in comparison with MIROC5 and observations. The 200-year-long data of the preindustrial control simulation by MIROC5 are used. The observations and reanalysis datasets used in the comparison are listed in Table 1.

Table 1 Summary of observation and reanalysis datasets used as references in the present paper.

Here, the model climatology in the preindustrial simulations is compared with observations in recent decades. Because observations are obtained concurrently with the progress of global warming due to increasing anthropogenic radiative forcing, the model climate under preindustrial conditions may not be adequate for use when making comparisons with recent observations. However, the root mean squared (RMS) errors of typical variables (e.g., the global mean SAT) in the climate models with respect to observations are much larger than the RMS differences between the model climatology in the preindustrial simulation and those in the last 30-year-long period in the historical simulations. Therefore, the differences between the time periods for which the climatology is defined are not a significant concern in comparisons among the climate models and observations.

3.1 Climatology

3.1.1 Atmosphere and land surface

First, model systematic biases in radiation at the TOA are evaluated because they reflect model deficiencies in cloud radiative processes that contribute to a large degree of uncertainty in climate modeling. Figure 4 shows annual mean biases in radiative fluxes at the TOA in MIROC6 and MIROC5 with respect to the recent Clouds and the Earth's Radiant Energy System (CERES) estimate (Loeb et al., 2009; the data are available at https://ceres.larc.nasa.gov/, last access: 12 June 2018). At the top right of each panel, a global mean (GM) value and a root mean squared error (RMSE) with respect to observations are written. In the present paper, RMSE is computed without model and observed global mean quantities unless otherwise noted.

Figure 4 Annual mean TOA radiative fluxes in MIROC6 (a, c, e) and MIROC5 (b, d, f). Upward is defined as positive. Net shortwave and longwave radiative fluxes and the sum of the two fluxes are denoted as OSR, OLR, and NET, respectively. Colors indicate errors with respect to observations (CERES) and contours denote values in each model. Note that a different color scale is used for longwave radiation. The global mean values and root mean squared errors are indicated by GM and RMSE, respectively. In the present paper, RMSE is computed without model and observed global mean quantities unless otherwise noted.

Persistent overestimates of net shortwave radiative flux and the sum of net shortwave and net longwave fluxes over low-latitude oceans in MIROC5 are significantly reduced in MIROC6. Hereafter, net shortwave radiation, net longwave radiation, and their sum are denoted as OSR, OLR, and NET, respectively, for simplicity. As described in Ogura et al. (2017), since parameter tuning cannot eliminate the abovementioned excess upward radiation, it is suggested that implementing a shallow convective parameterization is required in order to reduce the biases. Figure 5 shows annual mean moistening rates associated with deep and shallow convection at the 850 hPa pressure level in MIROC6. Moistening due to shallow convection occurs mainly over the low-latitude oceans, especially the eastern subtropical Pacific and the western Atlantic and Indian oceans. These active regions of shallow convection occur separately from regions with active deep convection in the western tropical Pacific and the ITCZ. The clear separation of the two convection types is consistent with satellite-based observations (Williams and Tselioudis, 2007). Owing to the shallow convective process that mixes the humid air in the planetary boundary layer with the dry air in the free troposphere, low-level cloud cover over the low-latitude oceans is better represented in MIROC6 than in MIROC5. Figure 6 shows annual mean biases in cloud covers with respect to the International Satellite Cloud Climatology Project (ISCCP; Rossow et al., 1996; Zhang et al., 2004; the data are available at https://isccp.giss.nasa.gov/, last access: 26 February 2018). An overestimate of low-level cloud cover over the low-latitude oceans in MIROC5 (Fig. 6b) is apparently reduced in MIROC6 (Fig. 6a), which results in smaller NET and OSR biases (Fig. 4). RMS error in low-level cloud cover in MIROC6 is 9 % lower than that in MIROC5.

Figure 5 Annual mean moistening rate associated with (a) deep convection and (b) shallow convection in MIROC6 at the 850 hPa pressure level.

Figure 6 Same as Fig. 4, but for cloud cover in MIROC6 (a, c, e) and MIROC5 (b, d, f). Low-, middle-, and high-level cloud cover is aligned from the top to the bottom. The tops for low-, middle-, and high-level clouds are defined to exist below the 680 hPa, between the 680 and 440 hPa, and above the 440 hPa pressure levels, respectively. The unit is nondimensional. ISCCP climatology is used as observations.

OSR in the midlatitudes is also better represented in MIROC6 than in MIROC5. Zonally distributed downward OSR bias in MIROC5 is reduced or becomes a relatively small upward bias in MIROC6 (Fig. 4c, d). This difference in the OSR bias is commonly found in both hemispheres. Cloud cover at middle and high levels is larger in MIROC6 over the subarctic North Pacific, North Atlantic, and the Southern Ocean (Fig. 6c–f), while low-level cloud cover over the same regions is smaller in MIROC6 than in MIROC5 over the same regions (Fig. 6a, b). The smaller low-level cloud cover in MIROC6 is inconsistent with the larger upward OSR bias in MIROC6. The wintertime midlatitude westerlies are stronger and are located more poleward in MIROC6 than in MIROC5. Correspondingly, the activity of sub-weekly disturbances in the midlatitudes is strengthened in MIROC6 (details are described later). These differences in the midlatitude atmospheric circulation between MIROC6 and MIROC5 lead to an enhanced poleward moist air transport from the subtropics to the subarctic region, which could result in an increase in the mid- and high-level cloud cover in MIROC6, as reported in previous modeling studies (e.g., Bodas-Salcedo et al., 2012; Williams et al., 2013). Consequently, the downward OSR bias in the midlatitudes is smaller in MIROC6 than in MIROC5. In polar regions, both biases in OSR and NET remain the same as in MIROC5.

Systematic bias in the outgoing longwave radiative flux (hereafter OLR) is worse in MIROC6 than in MIROC5 because MIROC6 tends to underestimate OLR over almost the entire global domain, except for Antarctica (Fig. 4e, f). The global mean of the high-level cloud cover in MIROC6 is larger than in MIROC5 by 0.04 (Fig. 6e, f), which is consistent with the smaller OLR in MIROC6. The increased moisture transport due to the strengthening of the westerlies and sub-weekly disturbances can partly explain the increase in the midlatitude high-level clouds in MIROC6, but high-level cloud cover is also larger in the low latitudes. Hirota et al. (2018) reported that moistening of the free troposphere due to shallow convection creates favorable conditions for atmospheric instabilities that lead to the resultant activation of deep convection at the low latitudes. Such processes may contribute to the inferior representation of OLR in MIROC6.

Next, we will discuss the global budget of the radiative fluxes and the RMS errors between models and observations. Note that only deviations from the global means are considered when calculating RMS errors. As shown in the upper right of Fig. 4a and b, the global mean (RMS errors) NETs are −1.11 (12.7) W m⁻² in MIROC6 and −0.98 (15.9) W m⁻² in MIROC5, respectively, and these values are consistent with the observed value of −0.81 W m⁻² (CERES; Loeb et al., 2009). However, the observed value is estimated in the present-day condition. Ideally, the model value in the preindustrial condition should be 0 W m⁻² and is in the marginally acceptable range. If NET is divided into OSR and OLR, so-called error compensation becomes apparent. The global means of OSR (OLR) are −231.3 (230.2) W m⁻² in MIROC6 and −237.6 (236.6) W m⁻² in MIROC5 (Fig. 4c–f). The observed global means of OSR and OLR are −240.5 and 239.7 W m⁻². Biases in the global mean OSR (OLR) with respect to observations are 9.2 (−9.5) W m⁻² in MIROC6 and 2.9 (−3.1) W m⁻² in MIROC5. Thus, the global mean OSR and OLR in MIROC6 are worse than those in MIROC5. Further division of OSR and OLR into cloud radiative forcing and clear-sky shortwave (longwave) radiative components shows that shortwave cloud radiative forcing is dominant on the biases in radiative fluxes. The biases in the global mean shortwave (longwave) cloud radiative forcing with respect to observations are 12.0 (6.7) W m⁻² in MIROC6 and −4.0 (−0.2) W m⁻² in MIROC5.

The global radiation budget in MIROC6 is inferior to that in MIROC5, while the reproducibility of the climatic means of typical model variables, other than radiative fluxes, and internal variations are better simulated in MIROC6 (details are shown later). As described in Sect. 2.5, the intensive tuning by perturbing model parameters is done by focusing on the reproducibility of climatic means, internal variations, and radiative forcing due to anthropogenic aerosols. During this procedure, the global radiation budget is traded off. On the other hand, RMS errors in NET, OSR, and OLR are 12.7, 16.2, and 6.3 W m⁻² in MIROC6 and 15.9, 18.9, and 6.8 W m⁻² in MIROC5, respectively, thereby indicating that the errors in MIROC6 have been reduced by 7 % to 20 %. This is also the case for shortwave and longwave cloud radiative forcings, for which the corresponding errors have been reduced by 17 % and 13 %, respectively. Taken together, these results show that the spatial patterns of the radiative fluxes are better simulated in MIROC6 than in MIROC5.

The improvement in spatial radiation patterns, especially in low-latitude OSR, is explained primarily by the implementation of shallow convective processes, which results in a moister free troposphere in MIROC6 than in MIROC5. Figure 7a and b show zonal mean biases in annual mean specific humidity with respect to the European Centre for Medium-Range Weather Forecasts Interim Reanalysis (ERA-I; Dee et al., 2011; the data are available at https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim, last access: 1 June 2019). Dry bias at 30^∘ S–30^∘ N, which occurs persistently in MIROC5, are largely reduced in MIROC6 owing to vertical mixing at the interface of the planetary boundary layer and the free troposphere. On the other hand, moist bias below the 600 hPa pressure level in the midlatitudes is somewhat worse in MIROC6 than in MIROC5. Shallow convection also contributes to the improvement of precipitation in the low latitudes. Figure 8 shows global maps for climatological precipitation in boreal winter (December–February) and summer (June–August). The second version of the Global Precipitation Climatology Project (GPCP; the data are available at https://precip.gsfc.nasa.gov/, last access: 21 December 2015) Monthly Precipitation Analysis (Adler et al., 2003) is used for the observations. While MIROC5 suffers from an underestimate of summertime precipitation over the western tropical Pacific, the underestimate is largely reduced in MIROC6 (Fig. 8d, f). The increase in precipitation is associated with deep convection because the moister free troposphere in MIROC6 is more favorable for the occurrence of deep convection (Hirota et al., 2018). On the other hand, the model representation of precipitation in MIROC6 is not necessarily alleviated other than the western tropical Pacific. For example, the overestimate of wintertime precipitation over the Indian Ocean and the midlatitude North Pacific is worse in MIROC6 than in MIROC5.

Figure 7 Annual and zonal mean specific humidity (a, b), temperature (c, d), and zonal wind (e, f) in MIROC6 (a, c, e) and MIROC5 (b, d, f). Colors indicate errors with respect to observations (ERA-I) and contours denote values in each model.

Figure 8 Precipitation in boreal winter (December–February; a, c, e) and summer (June–August; b, d, f) in observations (a, b; GPCP), MIROC6 (c, d), and MIROC5 (e, f). Areas with precipitation less than 3 mm d⁻¹ are not colored.

Zonal mean biases in annual mean air temperature and zonal wind velocity are also better represented in MIROC6 than in MIROC5 (Figs. 7c-f). The upper stratospheric warm bias at 50^∘ S–50^∘ N in MIROC5 is significantly reduced in MIROC6. The model top of MIROC6 is located at the 0.004 hPa pressure level and there are 42 vertical layers above the 50 hPa pressure level, while the model top of MIROC5 is placed at the 3 hPa pressure level. As a result, there are significant differences in stratospheric circulation between the models. As shown in the annual mean mass streamfunction calculated using zonal mean meridional winds (Fig. 9), an upward wind continuing from the low-latitude troposphere to the stratosphere is stronger in MIROC6 than in MIROC5. An increased upward advection of the temperature minimum around the tropopause at 30^∘ S–30^∘ N may lead to a reduction of warm temperature bias in the stratosphere, which is significant in MIROC5. Correspondingly, the stratospheric westerly bias at low latitudes of MIROC5 is also considerably alleviated in MIROC6. Note that the atmospheric O₃ concentration data used in MIROC5 are different from those in MIROC6, and the concentration in the stratosphere is higher than the data used in MIROC6. About 25 % of the abovementioned reduction in the stratospheric warm biases is explained by the smaller absorption of shortwave radiation by O₃. Note that the zonal mean temperature bias in Fig. 7c is smaller when the climatological mean temperature from 1980 to 2009 in a historical simulation is evaluated against observations because of the known stratospheric cooling with increased greenhouse gases and reduced O₃ concentrations.

Figure 9 Annual mean mass streamfunctions in (a) MIROC6 and (b) MIROC5. Contour interval is 0.3(0.025)×10¹⁰ kg s⁻¹ below (above) the 100 hPa pressure level. Negative values are denoted by dashed contours, and the horizontal dashed lines indicate the 100 hPa pressure level.

The zonal means of the air temperature and zonal wind in MIROC6 are also better simulated in the middle and high latitudes. A pair of easterly and westerly biases in MIROC5, which is in the troposphere of the Northern Hemisphere, is associated with a weaker midlatitude westerly jet and its southward shift with respect to observations. The pair of biases is reduced in MIROC6, thereby suggesting that a strengthening and northward shift of the westerly jet occur in MIROC6. Indeed, as shown in Fig. 10, the meridional contrast of high and low biases at the 500 hPa pressure level (Z500) along the wintertime westerly jet is weaker in MIROC6 than in MIROC5. The latitudes with the maximal meridional gradient of Z500 are located further northward in MIROC6 than in MIROC5, especially over the North Atlantic. Correspondingly, wintertime storm track activity (STA), which is defined as an 8 d high-pass-filtered eddy meridional temperature flux at the 850 hPa pressure level, is stronger over the North Pacific and Atlantic in MIROC6 than in MIROC5 (see Fig. 11) and is accompanied by an associated increase in precipitation, especially in the North Pacific (Fig. 8c, e). In the stratosphere above the 10 hPa pressure level, the polar night jet is reasonably captured in MIROC6, although the westerly is somewhat overestimated at 30–60^∘ N. Also, in the Southern Hemisphere, representation of the tropospheric westerly and the polar night jets is better in MIROC6 than in MIROC5, and the easterly bias centered at 60^∘ S in the troposphere is clearly reduced in MIROC6. Although causality is unclear, the warm air temperature bias above the tropopause to the south of 60^∘ S is smaller in MIROC6 than in MIROC5.

Figure 10 Same as Fig. 4, but for the wintertime 500 hPa pressure level in MIROC6 (a, c) and MIROC5 (b, d). Maps for boreal and austral winter are shown in (a, b) and (c, d), respectively. ERA-I is used as observations.

Figure 11 Wintertime storm track activity (STA) in observations (a, d), MIROC6 (b, e), and MIROC5 (c, f). STA is defined as 8 d high-pass-filtered eddy meridional temperature flux at the 850 hPa pressure level. Maps for boreal and austral winter are shown in (a–c) and (d–f), respectively. ERA-I is used as observations.

The enhanced wintertime STA in MIROC6 leads to a strengthening of the Ferrel circulation in the Northern Hemisphere and a broadening of its meridional width. As shown in Fig. 9, the northern edge of the Ferrel cell is located further northward in MIROC6 than in MIROC5. Because the Ferrel cell is a thermally indirect circulation driven primarily by eddy temperature and momentum fluxes, the stronger STA in MIROC6 possibly causes the Ferrel cell differences between the two models. Associated with the northward extension of the Ferrel cell, the upward wind between the Ferrel cell and the polar cell centered at 65^∘ N is stronger in MIROC6 than in MIROC5 and the meridional width of the polar cell is smaller. Also, in the Southern Hemisphere, the upward wind around 60^∘ S at the southern edge of the Ferrel cell is stronger in MIROC6 than in MIROC5. Correspondingly, high sea level pressure (SLP) biases in the polar region in MIROC5 are significantly reduced in MIROC6 (figures are omitted) and RMS errors with respect to observations (ERA-I) are decreased by 30 %. Meanwhile, in the stratosphere, anticlockwise (clockwise) circulation to the north (south) of 50^∘ N (S) is stronger and extends further upward in MIROC6 than in MIROC5. This circulation seems to continue from the troposphere into the stratosphere, thereby implying that more active troposphere–stratosphere interactions associated with wave coupling exist in MIROC6. Further details will be described later, focusing on the occurrence of sudden stratospheric warmings.

Parameterizations of SSNOWD (Liston, 2004; Nitta et al., 2014) and a wetland due to snowmelt water have been newly implemented into MIROC6 (Nitta et al., 2017). In comparison of MIROC6 with MIROC5, it can be seen that the former parameterization brings about significant improvement in the Northern Hemisphere snow cover fractions from the early to the late winter (Fig. 12). Compared with observations of the Northern Hemisphere EASE-Grid 2.0 (Brodzik and Armstrong, 2013; the data are available at https://nsidc.org/data/ease/, last access: 1 January 2013), the distribution of the snow cover fractions is more realistic in MIROC6 than MIROC5, especially where and when the snow water equivalent is relatively small (e.g., middle and high latitudes in November, over Siberia in February). Note that no clear improvement is found in May. This is because the newly implemented SSNOWD represents hysteresis in the relationship between snow water equivalent and snow cover fraction in both the accumulation and ablation seasons. MIROC6 underestimates the snow cover fraction in the partially snow-covered regions and overestimates it on the Tibetan Plateau and in some parts of China. We note that meteorological (e.g., precipitation or temperature) phenomena might affect these biases, but further investigation will be necessary to identify their causes. Nevertheless, in spite of those discrepancies, it can be said that the seasonal changes in the snow cover fraction are better simulated in MIROC6 than in MIROC5 (Fig. 12j).

Figure 12 Snow cover fractions for observations (a, d, g), MIROC6 (b, e, h), and MIROC5 (c, f, j). Maps in November, February, and May are aligned from the left to the right. The unit is nondimensional. Areas where snow cover fractions are less than 0.01 are masked. “Ave” and “corr” in the panels indicate spatial averages and correlation coefficients between observations and models over the land surface in the Northern Hemisphere, respectively. Time series in (j) shows the temporal rate of change of the monthly spatial averages. The snow cover dataset of the Northern Hemisphere EASE-Grid 2.0 is used as observations.

3.1.2 Ocean

Next, we evaluate the climatological fields of the ocean hydrographic structure, meridional overturning circulation (MOC), and sea ice distribution. The zonal mean potential temperature and salinity are displayed in Figs. 13 and 14, respectively. Both MIROC6 and MIROC5 capture the general features of the observed climatological hydrography (ProjD; Ishii et al., 2003). However, the potential temperatures in the deep and bottom layers to the south of 60^∘ S in the two models are warmer than observations because of insufficient formation and sinking of cold and dense water due to intense surface cooling around Antarctica (Figs. 13a–c and 14a–c). Such a warm temperature bias associated with deepwater formation is also found at the northern high latitudes of the Atlantic sector (Figs. 13a–c). By horizontal advection of the warm temperature biases associated with the Pacific and Atlantic MOC, the model temperatures in deep layers apart from polar regions are also warmer than in observations. The warm potential temperature bias in the deep layer is worse in MIROC6 than in MIROC5 in both the Atlantic and Pacific sectors, and the warm bias influences the subsurface and the intermediate layers above the 3000 m depth, which might be attributed to the excess ocean heat uptake and longer integration time in MIROC6 than in MIROC5 (the spin-up duration of MIROC6 is 2000 years and that of MIROC5 is about 1000 years). Also, the low-salinity bias below the 2000 m depth is worse in MIROC6 than in MIROC5, especially in the Pacific sector (Fig. 14e, f). This worsening can be explained by the excess supply of the freshwater in the Southern Ocean and weaker northward intrusion of the less saline water in MIROC6.

Figure 13 Annual mean potential temperature (a, b, c; unit is ^∘C) and salinity (d, e, f; psu) in the Atlantic sector for observations (a, d), MIROC6 (b, e), and MIROC5 (c, f). Colors indicate errors with respect to observations (ProjD) and contours denote model values in (b, e) and (c, f).

Figure 14 Same as Fig. 13, but for the Pacific sector.

In the Arctic Ocean, the halocline above the upper 500 m of depth is sharper and more realistic in MIROC6 than in MIROC5 and the high-salinity bias below the 500 m depth in MIROC5 is alleviated in MIROC6 (Fig. 13e, f) because, as described in Sect. 2.3, there are many more vertical levels in the surface and subsurface layers of MIROC6. In addition, vertical diffusivity in the Arctic Ocean is set to smaller values in MIROC6 than in MIROC5, and the turbulent kinetic energy input induced by surface wave breaking, as a function of the sea ice concentration in each grid cell, is reduced in MIROC6, as shown in Komuro (2014). In the North Pacific, the southward intrusion of North Pacific Intermediate Water (NPIW) around the 1000 m depth retreats northward in MIROC6. Strong tide-induced vertical mixing of seawater is observed along the Kuril Islands (e.g., Katsumata et al., 2004). The locally enhanced tide-induced mixing is known to reinforce the southward intrusion of the Oyashio and associated water mass transport from the subarctic to subtropical North Pacific and to feed the salinity minimum of NPIW (Nakamura et al., 2004; Tatebe and Yasuda, 2004). Hence, NPIW reproducibility is better in MIROC5, in which enhanced tidal mixing is considered, than in MIROC6. Because we encountered significant uncertainty in implementing the tidal mixing, we decided to stop implementing it in the development phase of MIROC6 at the expense of NPIW reproducibility.

The annual mean potential temperature and zonal currents along the Equator in MIROC6 are better simulated in MIROC6 than in MIROC5 (Fig. 15). Relatively cold water below the equatorial thermocline is upwelled in MIROC6, especially in the eastern tropical Pacific, which leads to a strengthening of the vertical temperature gradient across the thermocline. The eastward speed of the Equatorial Undercurrent in MIROC6 is over 80 cm s⁻¹, and is closer to the products of Simple Ocean Data Assimilation (SODA; Carton and Giese, 2008; the data are available at http://www.soda.umd.edu/, last access: 15 February 2019) than in MIROC5. These improvements are mainly attributed to the higher vertical resolution of MIROC6 in the surface and subsurface layers. However, the thermocline depths in the western tropical Pacific are still larger in the models than in observations and are attributed to the stronger trade winds in the models. When both MIROC6 and MIROC5 are executed as stand-alone AGCMs with the prescribed SST obtained from observations, the an overestimate of the equatorial trade winds also appears due to overestimate of the upward winds over the maritime continent associated with deep cumulus convection and the resultant strengthening of the Walker circulation over the equatorial Pacific. Better parameterizing deep cumulus convection in the models would be required for a better representation of the equatorial trade winds and thus oceanic states.

Figure 15 Annual mean climatology of potential temperature (^∘C; colors) and zonal current speed (cm s⁻¹; contours) along the Equator (1^∘ S–1^∘ N) in (a) observations (ProjD and SODA), (b) MIROC6, and (c) MIROC5.

Figure 16 displays annual mean Atlantic and Pacific MOC. In the Atlantic, two deep circulation cells associated with North Atlantic Deep Water (NADW; upper cell) and Antarctic Bottom Water (AABW; lower cell) are found in both of the models. NADW transport across 26.5^∘ N is 17.2 (17.6) Sv (1 Sv =10⁶ m³ s⁻¹) in MIROC6 (MIROC5). These values are consistent with the observational estimate of 17.2 Sv (McCarthy et al., 2015). RMS amplitudes of NADW transport are about 0.9 Sv in MIROC6 and 1.1 Sv in MIROC5 on longer than interannual timescales. These are smaller than the observed amplitude of 1.6 Sv in 2005–2014. Because observations include the weakening trend of the Atlantic MOC due to global warming, they can be larger than the model variability under preindustrial conditions. In the Pacific Ocean, both the models have the deep circulation associated with Circumpolar Deep Water (CDW), but the northward transport of CDW across 10^∘ S is 8.6 Sv in MIROC6, which is slightly larger than 7.5 Sv in MIROC5. Although these model values are somewhat smaller than observations, they are within the uncertainty range of observations (Talley et al., 2003; Kawabe and Fujio, 2010).

Figure 16 Annual mean meridional overturning circulation in the Atlantic (a, b) and the Indo-Pacific sectors (c, d) in MIROC6 (a, c) and MIROC5 (b, d). The unit is the sverdrup (Sv; ≡10⁶ m³ s⁻¹).

Northern Hemisphere sea ice concentrations are shown in Fig. 17. Here, it can be seen that both the March and September sea ice distributions in MIROC6 resemble the satellite-based observation (SSM/I; Cavarieli et al., 1991; the data are available at https://nsidc.org/, last access: 29 April 2019). In general, the spatial patterns of the models resemble the observations. Sea ice areas in March (September) are 12.4 (6.1), 13.0 (6.9), and 14.9 (5.7) million square kilometers in MIROC6, MIROC5, and observations, respectively. The model estimates are smaller (larger) in March (September) than in observations. The underestimate in March is still found in MIROC6 and is attributed to the underestimate of sea ice area in the Sea of Okhotsk and the Gulf of St. Lawrence, even though the sea ice area in the former region is better simulated in MIROC6 than in MIROC5. Meanwhile, the eastward retreat of the sea ice in the Barents Sea is better represented in MIROC6 than in MIROC5. The overestimates in September in the models are due to the model climatology being defined under preindustrial conditions, while observations are taken in present-day conditions of 1980–2009 when a rapid decreasing trend of summertime sea ice area (including a few events of drastic decreases) is ongoing (e.g., Comiso et al., 2008). Note that the model September sea ice area in 1980–2009 from historical simulations is smaller than the observations, and the sea ice area does not show a drastic year-to-year sea ice decrease with comparable amplitude with observations. The underestimate of the mean September sea ice area in MIROC6 might be attributed to slightly more rapid warming of the Arctic climate in MIROC6 than in observations. On the other hand, the modeled sea ice areas in the Southern Ocean are unrealistically smaller than in observations. Southern Hemisphere sea ice areas in March (September) are 0.1 (3.4), 0.2 (5.2), and 5.0 (18.4) million square kilometers in MIROC6, MIROC5, and observations, respectively. Since there are no significant differences between the two models, the spatial maps for the sea ice area in the Southern Hemisphere are omitted.

Figure 17 Northern Hemisphere sea ice concentrations in March (a, b, c) and September (d, e, f) for observations (a, d), MIROC6 b, e), and MIROC5 (c, f). The unit is nondimensional. Satellite-based sea ice concentration data of the SSM/I are used as observations.

Figure 18 shows the global maps of annual mean sea level height relative to the geoid. The absolute dynamic height values provided by Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO; Rio et al., 2014) data are used as observed sea level height (the data are available at https://www.aviso.altimetry.fr/en/home.html, last access: 20 June 2019). Overall, oceanic gyre structures in the two models are consistent with observations. Although representation of the gyres in MIROC6 remains generally the same as in MIROC5, there are a few improvements in the North Pacific and the North Atlantic. The midlatitude westerly in MIROC6 is stronger and is shifted further northward than in MIROC5 (Fig. 10), which results in the strengthening of the subtropical gyres, northward shifts of the western boundary currents, and their extensions. In particular, the current speed of the Gulf Stream and the North Atlantic Current is faster in MIROC6 than in MIROC5, and the contours emanating from the North Atlantic reach the Barents Sea in MIROC6. A corresponding increase in warm water transport from the North Atlantic to the Barents Sea leads to sea ice melting and an eastward retreat of the wintertime sea ice there in MIROC6 (Fig. 17a–c). An improvement in MIROC6 is also found in the Subtropical Countercurrent (STCC) in the North Pacific along 20^∘ N. As reported in Kubokawa and Inui (1999), the low-potential-vorticity water associated with a wintertime mixed layer deepening in the western boundary current region is transported southward in the subsurface layer, and it pushes up isopycnal surfaces around 25^∘ N. Thus, the eastward-flowing STCC is induced around 25^∘ N. Although both of the models show the wintertime mixed layer deepening, the ocean stratification along 160^∘ E is weaker in MIROC6 than in MIROC5 (not shown). This suggests that the isopycnal advection of low-potential-vorticity water in MIROC6 is more realistic than in MIROC5.

Figure 18 Annual mean sea level height relative to the geoid in (a) observations, (b) MIROC6, and (c) MIROC5. Contour interval is 20 cm. Negative values are denoted by dashed lines. Note that loading due to sea ice and accumulated snow on sea ice are removed from the model sea level height and that the global mean value is eliminated.

3.1.3 Discussions on model climatological biases

We have evaluated the simulated climatology in MIROC6 in comparison with MIROC5 and observations. The model climatology in MIROC6 shows certain improvements in simulating radiation, atmospheric and oceanic circulation, and the snow cover fractions in the Northern Hemisphere. In Fig. 19, we display the model biases in annual mean SAT and SST (Fig. 19) because these are typical variables that reflect errors in individual processes in the climate system. The global mean SAT (SST) is 15.2 (18.1) ^∘C in MIROC6, 14.6 (18.0) ^∘C in MIROC5, and 14.4 (18.1) ^∘C in observations. The modeled global mean SATs and SSTs are generally consistent with observations. However, since the observed (model) value is estimated in the present-day (preindustrial) condition, the model global mean SATs and SSTs are overestimated. Here, it should be noted that while the spatial patterns of the SAT and SST biases in MIROC6 resemble those in MIROC5, there are several improvements. For example, cold SAT bias in MIROC5 extending from the Barents Sea to Eurasia is significantly smaller in MIROC6, possibly owing to the increase in warm water transport by the North Atlantic Current and the resultant eastward retreat of the sea ice in the Barents Sea (Figs. 17 and 18). Warm SAT and SST biases along the west coast of North America are smaller in MIROC6 than in MIROC5. The reason is that an increase in southeastward Ekman transport in the eastern subarctic North Pacific due to the strengthening of the midlatitude westerly jet (Fig. 10) and the Aleutian low tend to cancel out the relatively warm water supply from the subtropics to the subarctic region by the surface geostrophic current. Although it is not clear from Fig. 19, the SAT and SST in the subtropical North Pacific around 20^∘ N are warmer by 2 K in MIROC6 than in MIROC5. Also in the Atlantic, the SAT in the western tropics is warmer in MIROC6. These warmer surface temperatures in MIROC6 indicate a reduction of the cold SAT and SST biases that can be alleviated by an increase in the downward OSR in MIROC6 due to the implementation of a shallow convective parameterization (Fig. 4) and by an increase in eastward transport of the warm pool temperature associated with the stronger STCC in MIROC6 (Fig. 18).

Figure 19 Same as Fig. 4, but for annual mean SAT (a, b) and SST (c, d). ERA-I for the SAT and the ProjD for the SST are used as observations.

On the other hand, the warm SAT and SST biases in the Southern Ocean and the warm SAT bias in the Middle East and the Mediterranean are worse in MIROC6 than in MIROC5. Consequently, the RMS error in SAT is larger in MIROC6 (2.4 K) than in MIROC5 (2.2 K). The former is essentially due to the underestimate of mid-level cloud cover, excess downward OSR, and the resultant underestimate of the sea ice in the Southern Ocean. Such a bias commonly occurs in many climate models and is normally attributed to errors in cloud radiative processes (e.g., Bodas-Salcedo et al., 2012; Williams et al., 2013). In addition, poor representations of mixed layer depths and open-ocean deep convection due to the lack of mesoscale processes in the Antarctic Circumpolar Current are causes of the warm bias (Olbers et al., 2004; Downes and Hogg, 2013). The latter warm bias, seen in the Middle East around the Mediterranean, can be explained by a tendency to underestimate the radiative forcing of aerosol–radiation interactions due to an underestimate of dust emissions from the Sahara in MIROC6 (not shown).

3.2 Internal climate variations

3.2.1 Madden–Julian oscillation and East Asian monsoon

In this section, we will evaluate the reproducibility of internal climate variations in MIROC6 in comparison with MIROC5 and observations, beginning with an examination of the equatorial waves in the atmosphere. Zonal wavenumber–frequency power spectra normalized by background spectra for the symmetric and antisymmetric components of OLR are calculated following Wheeler and Kiladis (1999) and are shown in Fig. 20. The daily mean OLR data derived from the Advanced Very High Resolution Radiometer (AVHRR) of the National Oceanic and Atmospheric Administration (NOAA) satellites (Liebmann and Smith, 1996; the data are available at https://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html, last access: 8 April 2019) are used for observational references. The signals corresponding to the Madden–Julian oscillation (MJO), equatorial Kelvin (EK), equatorial Rossby (ER), eastward inertia–gravity (n=1 EIG), and westward inertia–gravity (WIG) waves in the symmetric component and mixed Rossby–gravity (MRG) and eastward inertia–gravity (n=0 EIG) waves in the antisymmetric component stand out from the background spectra in observations. MIROC5 qualitatively reproduces these spectral maxima of the symmetric MJO, EK, and ER qualitatively, while the amplitudes of the MJO and the EK are underestimated. These underestimates are partially mitigated in MIROC6. The power summed over the eastward wavenumbers 1–3 and periods of 30–60 d corresponding to the MJO is 20 % larger in MIROC6 than in MIROC5. Furthermore, some additional analyses indicate that many aspects of the MJO, including its eastward propagation over the western tropical Pacific, are improved in MIROC6. Those improvements are primarily associated with the implementation of the shallow convective scheme that moistens the lower troposphere. The results of these additional analyses, along with some sensitivity experiments, are described in a separate paper (Hirota et al., 2018). The EIG and WIG in the symmetric component and the MRG and the EIG in the antisymmetric component are missing in both MIROC6 and MIROC5.

Figure 20 Zonal wavenumber–frequency power spectra of the (a–c) symmetric and (d–f) antisymmetric component of OLR divided by background power in (a, d) observations (NOAA OLR), (b, e) MIROC6, and (c, f) MIROC5. Dispersion curves of equatorial waves for the three equivalent depths of 12, 25, and 50 m are indicated by black lines. Signals corresponding to the westward and eastward inertia–gravity (WIG and EIG) waves, the equatorial Rossby (ER) waves, equatorial Kelvin waves, the mixed Rossby–gravity waves (MRG), and Madden–Julian oscillation (MJO) are labeled in (a). The unit of the vertical axes is cycles per day (cpd).

Figure 21 shows the June–August (JJA) climatology of precipitation and circulation in East Asia. As shown in observations (ERA-I; Fig. 21a), the East Asian summer monsoon (EASM) is characterized by the monsoon low over the warmer Eurasian continent and the subtropical high over the colder Pacific Ocean (e.g., Ninomiya and Akiyama, 1992). The southwesterly between these pressure systems transports moist air to the midlatitudes, forming a rainband called Baiu in Japanese. The general circulation pattern of the EASM and the rainband are well simulated in both MIROC6 and MIROC5. It should be noted that one of major deficiencies in MIROC5, the underestimate of the precipitation around the Philippines, has been largely alleviated in MIROC6. This improvement is, again, associated with the moistening of the lower troposphere by shallow convective processes. Interannual EASM variabilities are examined using an empirical orthogonal function (EOF) analysis of vorticity at the 850 hPa pressure level over [100–150^∘ E, 0–60^∘ N] following Kosaka and Nakamura (2010). The regressions of precipitation and 850hPa vorticity with respect to the time series of the first mode (EOF1) are shown in Fig. 21. In observations, precipitation and vorticity anomalies show a tripolar pattern with centers located around the Philippines, Japan, and the Sea of Okhotsk (Hirota and Takahashi, 2012). The anomalies around the Philippines and Japan correspond to the so-called Pacific–Japan pattern (Nitta et al., 1987). The southwest–northeast orientation of the wave-like anomalies is better simulated in MIROC6 than in MIROC5.

Figure 21 (a–c) Summertime (JJA) climatology of precipitation (colors; mm d⁻¹) and the 850 hPa horizontal wind (vector; m s⁻¹) for (a) observations (ERA-I), (b) MIROC6, and (c) MIROC5. (d–f) Anomalies of summertime precipitation (shading; mm d⁻¹) and the 850 hPa vorticity (contour; 10⁻⁶ s⁻¹) regressed to the time series of EOF1 of the 850 hPa vorticity over [100–150^∘ E, 0–60^∘ N] for (d) observations, (e) MIROC6, and (f) MIROC5.

Figure 22 shows the wintertime (December–February) climatology of circulation and the STA in East Asia. The East Asian winter monsoon (EAWM) is characterized by a northwesterly between the Siberian high and the Aleutian low in observations (ERA-I; e.g., Zhang et al., 1997). The monsoon northwesterly advects cold air to East Asia, enhancing the meridional temperature gradients and strengthening the subtropical jet around Japan. The jet's strength influences synoptic wave activities in the storm track. MIROC5 captures the circulation pattern but significantly underestimates the STA. The STA in MIROC6 is better simulated than in MIROC5, but it is still smaller than in observations. Interannual variability of the EAWM is also better represented in MIROC6 than in MIROC5. The dominant variability of the monsoon northwesterly is extracted as the EOF1 of the meridional wind at the 850 hPa pressure level over the region [30–60^∘ N, 120–150^∘ E]. In observations, the regressions with respect to the time series of the EOF1 show a stronger northwesterly accompanied by a suppressed STA, which is consistent with previous studies (Fig. 22d; e.g., Nakamura, 1992). This relationship between the circulation and the STA can be found in MIROC6 but not in MIROC5 (Fig. 22e, f). The explained variance of the EOF1 is 46.0 % in observations, 37.1 % in MIROC5, and 47.1 % in MIROC6, suggesting that the amplitude of this variability in MIROC6 is consistent with observations.

Figure 22 (a–c) Wintertime (DJF) climatology of STA (colors; K m s⁻¹), the 300 hPa zonal wind (contour; m s⁻¹), and the 300 hPa horizontal wind (vector; m s⁻¹) for (a) observations (ERA-I), (b) MIROC6, and (c) MIROC5. (d–f) As in (a–c), but for anomalies regressed onto the time series of the EOF1 of the 850 hPa meridional wind over [120–150^∘ E, 30–60^∘ N].

3.2.2 Stratospheric circulation

A few of the major changes in the model setting from MIROC5 to MIROC6 are higher vertical resolution and higher model-top altitude in MIROC6, namely the representation of the stratospheric circulation. Here, we examine the representation of the quasi-biennial oscillations (QBOs) in MIROC6. Figure 23 shows the time–height cross sections of the monthly mean, zonal-mean zonal wind over the Equator for observations (ERA-I) and MIROC6. In this figure, an obvious QBO with a mean period of approximately 22 months can be seen in MIROC6. The mean period is slightly shorter than that of ∼28 months in observations, and the simulated QBO period varies slightly from cycle to cycle. The maximum speed of the easterly at the 20 hPa pressure level is approximately −25 m s⁻¹ in MIROC6 and that of the westerly is 15 m s⁻¹. On the other hand, the observed maximum wind speeds are −35 m s⁻¹ for the easterly and 20 m s⁻¹ for the westerly. The simulated QBO has a somewhat weaker amplitude in MIROC6 than observations but the same east–west phase asymmetry. The QBO in MIROC6 shifts upward compared with that in observations, and the simulated amplitude is larger above the 5 hPa pressure level and smaller in the lower stratosphere. The simulated downward propagation of the westerly shear zones of zonal wind ($\partial \stackrel{\mathrm{‾}}{u}/\partial z>\mathrm{0}$, where z is the altitude) is faster than the downward propagation of easterly shear zones ($\partial \stackrel{\mathrm{‾}}{u}/\partial z)<\mathrm{0}$, which agrees with observations. The QBOs in MIROC6 are qualitatively similar to that represented in the MIROC ESM, which is an Earth system model with a similar vertical resolution that participated in CMIP5 (Watanabe et al., 2011). Note that nothing resembling a realistic QBO was simulated in the previous low-top version of MIROC5, which only has a few vertical layers in the stratosphere.

Figure 23 Time–height cross section of the monthly mean, zonal-mean zonal wind over the Equator for (a) observations (ERA-I) and (b) MIROC6. The contour intervals are 5 m s⁻¹. Dashed lines correspond to the altitude of the 70 hPa pressure level. Red and blue correspond to westerlies and easterlies, respectively.

Recently, Yoo and Son (2016) found that the observed MJO amplitude in the boreal winter is stronger than normal during the QBO easterly phase at the 50 hPa pressure level. They also showed that the QBO exerted greater influence on the MJO than did ENSO. Marshall et al. (2016) pointed out the improvement in forecast skill during the easterly phase of the QBO and indicated that the QBO could be a potential source of the MJO predictability. MIROC6 successfully simulates both the MJO and QBO in a way consistent with observations, as mentioned above, but correlations between the QBO and MJO are insignificant. One possible reason is the smaller amplitude of the simulated QBO in the lowermost stratosphere. The QBO contribution to tropical temperature variation at the 100 hPa pressure level is ∼0.1 K in MIROC6, which is much smaller than the observed value of ∼0.5 K (Randel et al., 2000). The simulated QBO has little effect on static stability and vertical wind shear in the tropical upper troposphere.

MIROC6 can also simulate sudden stratospheric warming (SSW), which is a typical intra-seasonal variability of the midlatitude stratosphere in the Northern Hemisphere. The standard deviation of monthly and zonal-mean zonal wind (colors) superimposed on a monthly climatology of zonal-mean zonal wind (black contours) in February is shown in Fig. 24a–c. There are two maxima of the standard deviations over the equatorial stratosphere and the middle- to high-latitude upper stratosphere in the Northern Hemisphere in observations (Fig. 24a), which correspond to QBO and polar vortex variability. This feature is well captured in MIROC6 (Fig. 24b), while there are too-small variations in MIROC5 where the stratosphere cannot be well resolved (Fig. 24c). The better representation of polar vortex variability in MIROC6 is closely associated with that of the SSW. As shown in Fig. 24, abrupt and short-lived warming events associated with SSW are detected in MIROC6, which are reproduced comparably to observations in terms of magnitude but are not detected in MIROC5. This is consistent with previous modeling studies that reported the importance of the well-resolved stratosphere for better simulation of stratospheric variability (e.g., Cagnazzo and Manzini, 2009; Charlton-Perez et al., 2013; Osprey et al., 2013). In December–January, however, MIROC6 still underestimates the frequency of SSW events, which is a common bias in other high-top climate models (e.g., Inatsu et al., 2007; Charlton-Perez et al., 2013; Osprey et al., 2013). It is conjectured that the less frequent SSW in December–January could be attributed to less frequent stationary wave breaking due to an overestimate of the climatological zonal wind speed of the polar night jet in MIROC6 (Fig. 24d and e).

Figure 24 Standard deviation of monthly and zonal-mean zonal wind (colors; unit is m s⁻¹) superimposed on monthly climatology of zonal-mean zonal wind (black contours; unit is m s⁻¹) in (a–c) February and (d–f) January for observations (ERA-I in 1979–2014; a, d, g), MIROC6 (b, e, h), and MIROC5 (c, f, i) during a 60-year period. In (g–i), the daily mean temperatures at the 10 hPa pressure level on the North Pole are plotted.

The inclusion of a well-resolved stratosphere in MIROC6 is also considered to be important for improvement in the representation of stratosphere–troposphere coupling. In order to evaluate this, we examine the time development of the northern annular modes (NAMs) associated with strongly weakened polar vortex events in the stratosphere. The NAM indices are defined by the first EOF mode of the zonal mean year-round daily geopotential height anomalies over the Northern Hemisphere and are computed separately at each pressure level (Baldwin and Thompson, 2009). The height anomalies are first filtered by a 10 d low-pass filter to remove transient eddies. Figure 25 shows the composite of the time development of the NAM index for weak polar vortex events. The events are determined by the dates on which the 10 hPa NAM index exceeded −3.0 standard deviations (Baldwin and Dunkerton, 2001). Note that the NAM index is multiplied by the square root of the eigenvalue in each level before the composite, that is, the composite having the geopotential height dimension. The weak polar vortex signal in the stratosphere propagates downward to the surface and persists for approximately 60 d in the lower stratosphere and upper troposphere. These observational features are well represented in MIROC6 (Fig. 25a, b). Although MIROC5 has also captured downward-propagating signals, its magnitude is approximately half in the stratosphere, and its persistency is weak in the lower stratosphere and upper troposphere. Therefore, these results strongly indicate that the inclusion of a well-resolved stratosphere in a model is important for representing not only stratospheric variability, but also stratosphere–troposphere coupling.

Figure 25 Composites of the time development of the zonal mean NAM index for stratospheric weak polar vortex events in (a) observations (ERA-I), (b) MIROC6, and (c) MIROC5. The indices have dimensions of geopotential height (m), and red denotes negative values. The color interval (contours) is 50 m (400 m). The number of events included in the composite is indicated above each panel.

3.2.3 El Niño–Southern Oscillation and Indian Ocean Dipole mode

Among the various internal climate variabilities on interannual timescales, ENSO is of great importance because it can influence climate not only in tropics, but also at the middle and high latitudes of both hemispheres through atmospheric teleconnections associated with wave propagations (e.g., Hoskins and Karoly, 1981; Alexander et al., 2002). Here, we describe the representation of ENSO and related teleconnection patterns. Figure 26 shows anomalies of SST, precipitation, the 500 hPa pressure height, and the equatorial ocean temperature regressed onto the Niño3 index, which is defined as the area average of the SST at [5^∘ S–5^∘ N, 150–90^∘ W]. ProjD and ERA-I in 1980–2009 are used as observations. Although the maximum of the SST anomalies in the tropical Pacific is shifted more westward than in observations, the ENSO-related SST anomalies simulated in both MIROC6 and MIROC5 are globally consistent with observations (Fig. 26a–c). Simulated positive precipitation anomalies in MIROC6 still overextend to the western Pacific (Fig. 26d–f). Meanwhile, dry anomalies over the maritime continent, the eastern equatorial Indian Ocean, and the SPCZ are better simulated in MIROC6 than in MIROC5. ENSO teleconnection patterns in Z500 (Fig. 26g–i) are also realistically simulated as seen in, for example, the Pacific–North American pattern (Wallace and Gutzler, 1981). Equatorial subsurface ocean temperature anomalies in MIROC6 are more confined within the thermocline than in MIROC5 (Fig. 26j–l), and the signals in MIROC6 are closer to observations. However, the subsurface signals in MIROC6 reside deeper than in observations. This is due to the difference in the climatological structure of the equatorial thermocline, which is attributed to the overestimate of the trade winds over the equatorial Pacific, as mentioned in Sect. 3.1.2.

Figure 26 Anomalies of SST (K), precipitation (mm d⁻¹), the 500 hPa pressure height (m), and the equatorial ocean temperature averaged over 5^∘ S–5^∘ N (K), which are regressed onto the Niño3 index. Monthly anomalies with respect to monthly climatology are used here. From the left to the right, the anomalies in observations (ProjD and ERA-I), MIROC6, and MIROC5 are aligned. Contours (j–l) denote annual mean climatological temperature with the 20 ^∘C isotherms thickened, and the contour interval is 2 ^∘C.

In addition to ENSO, the Indian Ocean Dipole (IOD) mode is recognized as a prominent interannual variability (Saji et al., 1999; Webster et al., 1999). Figure 27 shows anomalies of SST, 10 m wind, and precipitation regressed onto the autumn (September–November) dipole mode index (DMI), which is defined as the zonal difference of the anomalous SST averaged over [10^∘ S–10^∘ N, 50–70^∘ E] and that averaged over [10^∘ S–10^∘ N, 90–110^∘ E]. ProjD and ERA-I in 1980–2009 are used as observations. The observed positive IOD phase is characterized by a basin-wide zonal mode with positive (negative) SST anomalies in the western (eastern) Indian Ocean, and precipitation is increased (decreased) over the positive (negative) SST anomalies (Fig. 27a, d). The dipole SST pattern is better simulated in MIROC6 than in MIROC5 in which the eastern SST anomalies are located more southward than in observations (Fig. 27a–c). Correspondingly, a meridional dipole pattern in the precipitation of MIROC5 is alleviated, and MIROC6 shows a zonal dipole precipitation pattern, as in observations (Fig. 27d–f). Seasonal IOD phase locking to boreal autumn, which is assessed based on the RMS amplitude of the DMI, is also better simulated in MIROC6 than in MIROC5 (not shown). Seasonal shoaling of the eastern equatorial thermocline in the Indian Ocean is realistically simulated in MIROC6 during boreal summer to autumn. The shallower thermocline leads to the stronger thermocline feedback, which is evaluated based on the SST anomalies regressed onto the 20 ^∘C isotherm depth anomalies averaged over the eastern part of the IOD region. As displayed at the top of Fig. 27a–c, the thermocline feedback in MIROC6 is comparable to observations. This larger thermocline feedback in MIROC6 possibly leads to the abovementioned improvements in the IOD pattern. Note that the simulated surface wind anomalies are more realistic in MIROC6 than in MIROC5, although the magnitude of SST anomalies is overestimated in MIROC6. The overestimate of the SST anomalies may have arisen from an excessive response of the equatorial and coastal Ekman upwelling and down-welling to the wind changes, which are favorable in coarse-resolution ocean models.

Figure 27 Same as Fig. 26, but for anomalies of SST (colors), 10 m wind vectors (a, b, c), and precipitation (d, e, f) regressed onto the autumn DMI. The values of the regression slope between anomalies of the 20 ^∘C isotherm depth and the SST over the eastern IOD region, which indicates the thermocline feedback, are displayed on the top of (a, b, c).

3.2.4 Decadal-scale variations in the Pacific and Atlantic Ocean

On longer than interannual timescales, the PDO (Mantua et al., 1997) or the Interdecadal Pacific Oscillation (IPO; Power et al., 1999) is known to be a dominant climate mode that is detected in the SST and the SLP over the North Pacific. To examine simulated PDO patterns, monthly SST and wintertime (December–February) SLP anomalies are regressed onto the PDO index, defined as the first EOF mode of the North Pacific SST to the north of 20^∘ N, and shown in Fig. 28. In order to detect decadal-scale variation, the COBE-SST2–SLP2 data (Hirahara et al., 2014) from 1900 to 2013 are used as observations. Negative SST anomalies in the western and central North Pacific and positive SST anomalies in the eastern North Pacific are found in observations. These signals are also represented in both MIROC6 and MIROC5. The regression of SLP anomalies corresponding to the deepening of the Aleutian low are well simulated in the models over the subarctic North Pacific, and it can be seen that the amplitudes of the SLP anomalies are larger in MIROC6 than in MIROC5, which is closer to the observation. In the tropical Pacific, positive SST anomalies, which are among the more important driving processes of the PDO (e.g., Alexander et al., 2002), are seen in both the models and the observations. In MIROC5, the 5-year running means of the wintertime (November–March) North Pacific Index (NPI), defined as the SLP averaged over [30–65^∘ N, 160^∘ E–140^∘ W], are less sensitive to the Niño3 index (correlation coefficient $r=-\mathrm{0.37}$) than to the Niño4 index ($r=-\mathrm{0.64}$). Note that the Niño4 index is defined as the area average of the SST over [5^∘ S–5^∘ N, 160^∘ E–150^∘ W]. The distorted response of the extratropical atmosphere to the tropical SST variations works to unsuitably modify the extratropical ocean and plays a major role in limiting the decadal predictability of the PDO index in MIROC5 (Mochizuki et al., 2014). In contrast, those in MIROC6 are well correlated with the Niño3 index ($r=-\mathrm{0.61}$) in addition to the Niño4 index ($r=-\mathrm{0.62}$). Overestimates of the tropical signals of MIROC5 in the western tropical Pacific are also alleviated in MIROC6. The abovementioned PDO improvement and the linkage between the tropics and the midlatitude North Pacific imply potential for improved skills in initialized decadal climate predictions.

Figure 28 Same as Fig. 26, but for anomalies of monthly SST and wintertime SLP regressed onto the PDO index (see the text). COBE-SST2–SLP2 data in 1900–2013 are used as observations.

In the Atlantic Ocean, there is another decadal-scale variability, which is called the AMO (Schlesinger and Ramankutty, 1994). Figure 29 shows anomalies of SST and SLP regressed onto the AMO index, which is defined as the area average of the SST anomalies in the North Atlantic over [0–60^∘ N, 0–80^∘ W] with the global mean SST anomalies subtracted (Trenberth and Shea, 2006). As in the PDO, the centennial COBE-SST2–SLP2 data in 1900–2013 are used as observations. The observed AMO spatial pattern in its positive phase is characterized by positive SST anomalies in the off-Equator and the subarctic North Atlantic and by negative or weakly positive SST anomalies in the western subtropical North Atlantic (Fig. 29a). Corresponding to negative (positive) SLP anomalies over the subtropical (subarctic) North Atlantic, the midlatitude westerly jet is weaker in a positive AMO phase than in normal years. These spatial patterns in the SST and SLP are simulated in both MIROC6 and MIROC5. It is especially noteworthy that the positive SST anomalies at low latitudes have larger amplitudes in MIROC6 than in MIROC5, and they extend to the South Atlantic as in observations (Fig. 29b, c). On the other hand, the positive SST anomalies in the subarctic region are underestimated in MIROC6, which may be due to the smaller RMS amplitudes of NADW transport in MIROC6 (see Sect. 3.1).

Figure 29 Same as Fig. 26, but for anomalies of SST (colors) and SLP (contours; 0.2 hPa) regressed onto the AMO index (see the text). Negative values are denoted by dashed contours.

3.3 Climate sensitivity

Following the regression method by Gregory et al. (2004) and Gregory and Webb (2008), we conducted abrupt CO₂ quadrupling experiments with MIROC6 (Tatebe and Watanabe, 2018b) and MIROC5 in order to evaluate effective climate sensitivity (ECS), radiative forcing, and climate feedback. The CO₂ quadrupling experiments were initiated from the preindustrial control runs. Data from the first 150 years after the CO₂ increase were used for the analysis.

ECS, 2× CO₂ radiative forcing, and climate feedback for MIROC6 are estimated to be 2.6 K, 3.7 W m⁻², and −1.4 W m⁻² K⁻¹, respectively (Fig. 30a and Table 2). The ECS, radiative forcing, and climate feedback in MIROC6 are lower, higher, and negatively larger than those of the CMIP5 multi-model ensemble means, although these estimates for MIROC6 are within the ensemble spread of the multi-models (Andrews et al., 2012). The ECS of MIROC6 is almost the same as MIROC5 because the decrease in radiative forcing is counterbalanced by the positive increase in climate feedback, although the change in climate feedback is small and not statistically significant. The decrease in radiative forcing of MIROC6 relative to MIROC5 is evident in the longwave and shortwave cloud components (LCRE and SCRE in Fig. 30b and Table 3). On the other hand, the clear-sky shortwave component (SWclr) increases in MIROC6 relative to MIROC5, which partially cancels the differences between the two models. The positive increase in climate feedback is pronounced in the SCRE, which is partially offset by the decrease in the clear-sky longwave (LWclr) and SWclr (Fig. 30c and Table 3).

Table 2 Effective climate sensitivity (ECS), radiative forcing of CO₂ doubling, and climate feedback for MIROC6 and MIROC5. The MIROC6 result marked with an asterisk (^*) is different from MIROC5 at the 5 % level.

Table 3 Radiative forcing of CO₂ doubling and climate feedback for MIROC6 and MIROC5, evaluated with different components of TOA radiation as longwave clear sky (LWclr), shortwave clear sky (SWclr), longwave cloud radiative effect (LCRE), and shortwave cloud radiative effect (SCRE). The MIROC6 results marked with an asterisk (^*) are different from MIROC5 at the 5 % level.

Figure 30 (a) Global mean net radiative imbalance at the TOA plotted against the global mean SAT increase. Data from the first 150 years after the abrupt CO₂ quadrupling are used. (b) 2× CO₂ radiative forcing estimated by regressing four components of TOA radiation against the global mean SAT, following Gregory and Webb (2008). (c) Same as (b) but for climate feedback. In (b, c) LWclr (SWclr) and LCRE (SCRE) denote a clear-sky longwave (shortwave) component and a longwave (shortwave) cloud component, respectively. The arrows in (b) and (c) indicate that the results of MIROC6 are different from MIROC5 at the 5 % level.

We now focus on the SCRE of the radiative forcing and climate feedback, which show the largest differences between the two models, and compare the geographical distribution (Fig. 31). The distribution is calculated by regressing the changes in SCRE caused by the CO₂ increase at each latitude–longitude grid box against the change in the global mean SAT. There is a large difference in the geographical distribution between MIROC6 and MIROC5, with the former showing a more pronounced zonal contrast in the tropical Pacific than the latter. The changes in the global mean from MIROC5 to MIROC6 (Fig. 30b, c) are correlated with changes in the western tropical Pacific, showing more negative radiative forcing and more positive climate feedback, which are partially offset by the changes in the central tropical Pacific with opposite signs. The radiative forcing and climate feedback tend to show similar geographical patterns with opposite signs in each model.

Figure 31 Shortwave cloud component of (a, c) 2× CO₂ radiative forcing and (b, d) climate feedback in MIROC6 (a, b) and MIROC5 (c, d).

4 Summary and discussion

The sixth version of a climate model, MIROC6, was developed by a Japanese climate modeling community aiming to contribute to CMIP6 through a deeper understanding of a wide range of climate science issues and seasonal to decadal climate predictions and future climate projections. The model configurations and basic performance in the preindustrial control simulation have been described and evaluated in the present paper. Major changes from MIROC5, which was our official model for CMIP5, to MIROC6 are mainly done in the atmospheric component. These include the implementation of a parameterization of shallow convective processes, the higher model top, and vertical resolution in the stratosphere. The ocean and land surface components have been also updated in terms of the horizontal grid coordinate system and higher vertical resolution in the former, as well as parameterizations for sub-grid-scale snow distribution and wetlands due to snowmelt water in the latter. Overall, the model climatology and internal climate variability of MIROC6, which are assessed in comparison with observations, are better simulated than in MIROC5.

An overestimate of low-level cloud amounts at low latitudes, which can be partly attributed to insufficient representation of shallow convective processes, is significantly alleviated in MIROC6. The free atmosphere becomes wetter and the precipitation over the western tropical Pacific becomes larger in MIROC6 than in MIROC5, primarily due to vertical mixing of the humid air in the planetary boundary layer with the dry air in the free troposphere. Shallow convection also contributes to better propagation characteristics of intra-seasonal variability associated with MJO in MIROC6, as well as East Asian summer monsoon variability on interannual timescales. In addition, QBO, which is absent in MIROC5, appears in MIROC6 because of its better stratospheric resolution and non-orographic gravity wave drag parameterization.

Climatic mean and internal climate variability in the midlatitudes are also improved in MIROC6. Together with enhanced activity of sub-weekly disturbances, the tropospheric westerly jets in MIROC6 are shifted more poleward and are stronger than in MIROC5, especially in the Northern Hemisphere. Overestimates in zonal wind speed of the polar night jet are reduced in MIROC6. These advanced representations lead to tighter interactions between the troposphere and the stratosphere in MIROC6. SSW events, in the form of polar vortex destruction induced by upward momentum transfer from the troposphere to the stratosphere (e.g., Matsuno, 1971), are well captured in MIROC6. On interannual timescales, the improvement of the westerly jet results in better representations of the spatial wind pattern of the wintertime East Asian monsoon. Associated with changes in large-scale atmospheric circulation, the western boundary currents in the oceans, the Kuroshio–Oyashio current system, the Gulf Stream, and their extensions are better simulated in MIROC6. The increase in warm water transport from the subtropical North Atlantic to the Barents Sea seems to melt the sea ice in the Barents Sea and to alleviate the overestimate of the wintertime sea ice area that is seen in that region in MIROC5. Another improvement in MIROC6 is found in the climatological snow cover fractions in the early winter over the Northern Hemisphere continents. In the Southern Hemisphere, however, the underestimate of mid-level clouds and the corresponding warm SAT bias, the underestimate of sea ice area, and the overestimate of incoming shortwave radiation in the Southern Ocean, all of which are attributed to errors in cloud radiative and planetary boundary layer processes (e.g., Bodas-Salcedo et al., 2012; Williams et al., 2013), remain the same as in MIROC5.

Qualitatively, the linkage representations between the tropics and the midlatitudes associated with ENSO in MIROC6 are mostly the same as in MIROC5. Meanwhile, oceanic subsurface signals, which partly control ENSO characteristics, are more confined along the equatorial thermocline in MIROC6, which is consistent with observations. Regarding the PDO, tropical influence on the midlatitudes is more dominant in MIROC6 than in MIROC5, suggesting improvements in decadal-scale atmospheric teleconnections in MIROC6.

The above descriptions are mainly on the Pacific internal climate variabilities. Regarding the Indian Ocean, the zonal dipole structures in the SST and precipitation associated with the interannual variability, known as the IOD, are better simulated in MIROC6 than in MIROC5, which has a bias of a false meridional precipitation pattern. In the Atlantic, the multi-decadal variability, known as the AMO, is represented in both the models as roughly consistent with observations, but their reproducibility shows both drawbacks and advantages. Signals associated with AMO in the subarctic (tropical) region are underestimated in MIROC6 (MIROC5).

As a metric for climate change induced by atmospheric CO₂ increase, ECS is also estimated. Although the model configurations and performances are different between the models, the ECS is almost the same (2.6 K). However, looking at the geographical distributions of radiative forcing and climate feedback, the amplitudes of shortwave cloud components are much larger in MIROC6 than in MIROC5. Since the larger negative (positive) radiative forcing and positive (negative) climate feedback in the western (central) tropical Pacific cancel each other, global mean quantities in MIROC6 remain almost the same as in MIROC5. As a topic of future study, estimating radiative forcing and climate feedback with experiments like those in the Atmospheric Model Intercomparison Project in order to check the robustness of the present study would be desirable. Elucidating the impact of different geographical patterns of radiative forcing and climate feedback on the projected future climates would also be useful.

After conducting the preindustrial control simulation and evaluating the model reproducibility of the mean climate and the internal climate variability, ensemble historical simulations that were initiated from the preindustrial simulations were executed using the historical forcing data recommended by the CMIP6 protocol (Tatebe and Watanabe, 2018c). Figure 32 shows a time series of the global mean SAT anomalies with respect to the 1961–1990 mean. There are 30 (5) ensemble members in the MIROC6 (MIROC5) historical simulations. Note that the MIROC5 historical simulations are executed using the forcing datasets of the CMIP5 protocol. As shown in Fig. 32, the simulated SAT variations in both MIROC6 and MIROC5 follow observations (HadCRUTv4.4.0; Morice et al., 2012; the data are available at https://crudata.uea.ac.uk/cru/data/temperature/, last access: 5 June 2019) on a centennial timescale. The temperature increases from the nineteenth century to the early twenty-first century are about 0.72 K in MIROC6, 0.85 K in MIROC5, and 0.82 K in observations. Focusing on the period from the 1940s to the 1960s, the SAT variations seem to be better simulated in MIROC6 than in MIROC5, which can be due to both an update of the forcing datasets and the larger ensemble number in MIROC6. On the other hand, the warming trend during the first half of the twentieth century in the models is about half as large as in observations. Whether it can be attributed to internal climate variability (e.g., Thompson et al., 2014; Kosaka and Xie, 2016) or to an externally forced mode (e.g., Meehl et al., 2003; Nozawa et al., 2005) is still being debated. The so-called recent hiatus of global warming (Easterling and Wehner, 2009) in the first decade of the twenty-first century is not simulated in MIROC6 or MIROC5. The observed hiatus is considered to occur in association with a negative IPO phase as internal climate variation (e.g., Meehl et al., 2011; Watanabe et al., 2014). As external drivers of the hiatus, the increase in stratospheric water vapor and the weakening of solar activity are given as possible candidates (e.g., Solomon et al., 2010; Kaufmann et al., 2011). Failure to simulate the hiatus in the models could be attributed to uncertainties in the historical forcing datasets or cancellation of internal climate variations of the IPO by ensemble mean manipulation of the individual historical simulations.

Figure 32 Time series of the global mean SAT anomalies for observations (black), MIROC6 (red), and MIROC5 (blue). A 5-year running-mean filter is applied to the anomalies with respect to the 1961–1990 mean. Colors indicate spreads of ensemble experiments for each model (1 standard deviation).

As summarized above, the overall reproducibility of the mean climate and the internal variability in the latest version of our climate model, MIROC6, has progressed, as has the historical warming trend of the climate system. During the first trial of the preindustrial simulation conducted just after the model configuration was frozen, however, the model reproducibility was not as good as seen in MIROC5. As described in Sect. 2.5, we intensively tuned the model by perturbing parameters associated with cumulus and shallow convection and planetary boundary processes in particular. In addition, before starting the historical simulations, we estimated and tuned the radiative forcing due to aerosol–radiation and aerosol–cloud interactions by changing the parameters of cloud microphysics in order to ensure that the estimated radiative forcing would be closer to the best estimate of the IPCC AR5 (IPCC, 2013). Without this parameter tuning, the simulated warming trend after the 1960s was 70 % as large as seen in observations. This dependence of radiative forcing and reproducibility of the warming trend on cloud microphysics has also been reported in other climate models (Golaz et al., 2013). A recent comparison of cloud microphysical statistics between climate models and satellite-based observations has pointed out that “tuned” model parameters that were adjusted for adequate radiative forcing and realistic SAT changes do not necessarily ensure that cloud properties and rain–snow formation will be consistent with observations and implies the presence of error compensations in climate models (e.g., Suzuki et al., 2013; Michibata et al., 2016). Error compensation is also found in both global and regional aspects. As described in Sect. 3.1, the global TOA radiation imbalance in MIROC6 is about −1.1 W m⁻², which is in the acceptable range of observations. However, when the TOA imbalance is examined in parts, cloud radiative components in the model contain non-negligible biases with respect to satellite-based observations. Regarding error compensation in the oceanic processes, the modeled northward transport of CDW, which is within the uncertainty range of observations, is maintained by spurious open-ocean convection in the Southern Ocean that often appears in coarse-resolution ocean models for which oceanic mesoscale eddies and coastal bottom water formation cannot be represented (e.g., Olbers et al., 2004; Downes and Hogg, 2013).

Several key foci remain for ongoing model development efforts. These include process-oriented refinements of cloud microphysics and convective systems based on constraints from satellite data and feedbacks from cloud-resolving atmospheric models (e.g., Satoh et al., 2014), higher resolutions for representations of regional extremes, oceanic eddies and river floods, and parameterization of tide-induced microscale mixing of seawater. Improvement of computational efficiency, especially on massively parallel computing systems, is among the urgent issues for long-term and large-ensemble simulations. These improvements can contribute to a deeper understanding of the Earth's climate, reducing uncertainties in climate projections and predictions, and more precise evaluations of human influences on carbon–nitrogen cycles when applied to Earth system models.

Appendix A

Table A1 Summary of the updated configurations from MIROC5 to MIROC6.

NA: not available.