Improving Madden–Julian oscillation simulation in atmospheric general circulation models by coupling with a one-dimensional snow–ice–thermocline ocean model

A one-column, turbulent, and kinetic-energy-type ocean mixed-layer model (snow–ice–thermocline, SIT), when coupled with three atmospheric general circulation models (AGCMs), yields superior Madden–Julian oscillation (MJO) simulations. SIT is designed to have fine layers similar to those observed near the ocean surface; therefore, it can realistically simulate the diurnal warm layer and cool skin. This refined discretization of the near-surface ocean in SIT provides accurate sea surface temperature (SST) simulation, and thus facilitates realistic air–sea interaction. Coupling SIT with the European Centre/Hamburg Model version 5, the Community Atmosphere Model version 5, and the High-Resolution Atmospheric Model significantly improved MJO simulation in three coupled AGCMs compared to the AGCM driven by a prescribed SST. This study suggests two major improvements to the coupling process. First, during the preconditioning phase of MJO over the Maritime Continent (MC), the often underestimated surface latent heat bias in AGCMs can be corrected. Second, during the phase of strongest convection over the MC, the change in intraseasonal circulation in the meridional circulation enhancing near-surface moisture convergence is the dominant factor in the coupled simulations relative to the uncoupled experiments. The study results show that a fine vertical resolution near the surface, which better captures temperature variations in the upper few meters of the ocean, considerably improves different models with different configurations and physical parameterization schemes; this could be an essential factor for accurate MJO simulation.

1 Introduction

The Madden–Julian Oscillation (MJO) is the dominant pattern of atmospheric intraseasonal variability in the tropics (Madden and Julian, 1972; Zhang, 2005; Jiang et al., 2020). It has been reported that the MJO convection is often observed over sea surface temperature (SST) of greater than 28 ^∘C in the Indo-Pacific warm pool (Salby and Hendon, 1994). MJO is an eastward-propagating ocean–atmosphere and convection–circulation coupled phenomenon that lasts for 20–100 d. On these timescales, low-level moisture convergence, warm SST, and shallow upper-ocean mixed-layer depth precede the eastward propagation of organized deep convection by approximately 10 d, with opposite conditions following approximately 10 d afterwards (Krishnamurti et al., 1988; Hendon and Salby, 1994; Woolnough et al., 2000). Heat flux exchange between the atmosphere and ocean modulates the intraseasonal oscillation (Shinoda and Hendon, 1998). Studies have emphasized the importance of moisture and heat flux feedback in MJO (Sobel et al., 2008, 2010; DeMott et al., 2015). Besides, oceanic wave dynamics are suggested to be associated with MJO (e.g., zonal wind stress anomalies driven by the MJO-forced, eastward-propagating, oceanic equatorial Kelvin wave; Hendon et al., 1998; Webber et al., 2010), and the signals could extend as deep as 1500 m in the ocean (Matthews et al., 2007). Furthermore, the westward-propagating, oceanic, equatorial Rossby wave can initiate the next MJO in the Indian Ocean (Webber et al., 2010, 2012). Evidence of oceanic intraseasonal signals coupling with atmospheric signals was observed in terms of the sea level, surface heat flux, salinity, and temperature during field experiments and in situ monitoring (Oliver and Thompson, 2011; Drushka et al., 2012; Wang et al., 2013; Chi et al., 2014; Matthews et al., 2014; DeMott et al., 2015; Fu et al., 2015).

Recent modeling studies have demonstrated that most coupled models could improve MJO simulations, but the ocean, driven by the atmosphere, contributes indirectly by improving the mean state, heat flux, fresh water, and momentum. DeMott et al. (2016) estimated that direct SST-driven ocean feedback contributes to the MJO propagation of up to 10 % by a change in column moisture. A comparison of the direct and indirect effects of SST indicated that direct effects, such as SST-driven surface fluxes, tend to offset wind-driven fluxes (DeMott et al., 2015, 2016, 2019). The factor of indirect ocean feedback on the atmospheric physical process includes strong MJO convection, which can amplify the radiative feedback to MJO convections associated with large cloud systems (Del Genio and Chen, 2015). The SST gradients can drive the MJO low-level convergence (Hsu and Li, 2012; Li and Carbone, 2012) and destabilize lower tropospheric conditions to further enhance low-level convergence to the east of the MJO convergence (Wang and Xie, 1998; Marshall et al., 2008; Benedict and Randall, 2011; Fu et al., 2015). Many observational and model studies have reported that coupled feedback enhances the MJO with strong horizontal moisture advection, driven by sharp mean near-equatorial meridional moisture gradients (DeMott et al., 2015; Jiang et al., 2018; DeMott et al., 2019; Jiang et al., 2020). These findings suggest that high-frequency SST perturbations could improve moisture convergence efficiency and enhance MJO propagation through relatively smooth background moisture distribution.

Tseng et al. (2015) identified the key role of the upper-ocean warm layer in improving the MJO eastward propagation simulation using the European Centre/Hamburg Model version 5 (ECHAM5), coupled with the one-column ocean model named snow–ice–thermocline (SIT). Many observational (Drushka et al., 2012; Chi et al., 2014) and modeling (Klingaman and Woolnough, 2013; DeMott et al., 2019; Klingaman and Demott, 2020) studies have supported this hypothesis. However, coupling the SIT to only one atmospheric general circulation model (AGCM) may be insufficient to prove the effectiveness of the coupling. In this study, we coupled the SIT to three AGCMs: ECHAM5, the Community Atmosphere Model version 5 (CAM5), and the High-Resolution Atmospheric Model (HiRAM). We also discussed the coupling mechanism that leads to simulation improvement.

The remainder of the paper is organized as follows. In Sect. 2, we describe the models, experimental designs, and observational data. Sections 3 and 4 present the results and discussion, respectively.

2 Data, model experiments, and methodology

2.1 Observational, atmospheric, and oceanic data

Observational data used in this study include precipitation from the Global Precipitation Climatology Project version 1.3 (GPCP, 1^∘ resolution, 1997–2010; Adler et al., 2003), outgoing longwave radiation (OLR, 1^∘ resolution, 1997–2010; Liebmann, 1996), and daily SST (optimum interpolated SST, OISST, 0.25^∘ resolution, 1989–2010; Banzon et al., 2014.) from the National Oceanic Atmosphere Administration. The in situ ocean temperature profiles from 1989 to 2010 were obtained from the Tropical Ocean Global Atmosphere program (McPhaden et al., 2010).

Atmospheric variables were obtained from the European Centre for Medium-range Weather Forecast Reanalysis (ERA-Interim, Dee et al., 2011) from 1989 to 2010. The variables include zonal wind, meridional wind, temperature, specific humidity, sea level pressure, geopotential high, latent heat, sensible heat, and shortwave and longwave radiation. Oceanic temperature data from 1989 to 2010 were obtained from the NCEP Global Ocean Data Assimilation System (GODAS) (Behringer and Xue, 2004) provided by NOAA/OAR/ESRL PSL (Boulder, Colorado, USA; https://psl.noaa.gov/data/gridded/data.godas.html, last access: 11 August 2020).

2.2 Model experiments

In this study, we coupled the one-column ocean model SIT (Tu and Tsuang, 2005; Tsuang et al., 2009) to three AGCMs. SIT simulates variations in the SST and upper-ocean temperature, including the diurnally varying cool skin and warm layer in the upper few meters of the ocean and the turbulent kinetic energy (TKE; Gaspar et al., 1990) in the water column (Tsuang et al., 2001; Tu and Tsuang, 2005; Tu, 2006; Tsuang et al., 2009; Tseng et al., 2015; Lan et al., 2021a​​​​​​​). Cool skin is a very thin layer that has a direct contact with the atmosphere, and the warm layer is the warmer sea water immediately below the cool skin in the top few meters of the ocean. They fluctuate diurnally in response to atmospheric forcing. SIT with high vertical resolution realistically simulates the warm layer (within top 10 m) and cool skin (the top layer with 0.001 m thickness) and improves the simulation of upper-ocean temperature (Tu and Tsuang, 2005; Tsuang et al., 2009). The model has been verified at a tropical ocean site (Tu and Tsuang, 2005), in the South China Sea (Lan et al., 2010), and Caspian Sea (Tsuang et al., 2001). The melt and formation of snow and ice above a water column have also been introduced (Tsuang et al., 2001). The three AGCMs used in this study are as follows. ECHAM5, the fifth-generation AGCM developed at the Max Planck Institute for Meteorology (Roeckner et al., 2003, 2006), is a spectral model that employs the Nordeng cumulus convective scheme (Nordeng, 1994). We used a horizontal resolution of T63 (approximately 2^∘) with 31 vertical layers and a model top at 10 hPa (approximately 30 km). The second model is the NCAR Community Atmospheric Model version 5 (Hurrell et al., 2013) from the National Center for Atmospheric Research. We used a horizontal resolution of approximately 1.875^∘ latitude × 2.5^∘ longitude and 30 vertical layers with the Zhang–McFarlane deterministic convection scheme (Zhang and McFarlane, 1995) and the University of Washington shallow convection scheme (Park and Bretherton, 2009). HiRAM was developed based on the Geophysical Fluid Dynamical Laboratory global atmosphere and land model (AM2, Team et al., 2004; Zhao et al., 2009) with a few modifications (Chen et al., 2019). We used a horizontal resolution of 0.5^∘ latitude × 0.5^∘ longitude with 32 vertical levels. For the boundary layer and free atmospheric turbulence, the model adopted the 2.5-order parameterization of Mellor and Yamada (1982). Surface fluxes are computed based on the Monin–Obukhov similarity theory, given the atmospheric model's lowest level of wind, temperature, and moisture.

There are 42 vertical layers in SIT, with 12 layers in the upper 10 m: the surface, 0.05 mm, 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 m below the ocean surface. The fine resolution was designed to realistically simulate the upper-ocean warm layer, including a layer at 0.05 mm, reproducing the cool skin of the ocean surface. It is worth noting that coupling a high vertical resolution TKE ocean model with an AGCM is unconventional. To account for neglected horizontal processes, the model ocean was weakly nudged (with a 30 d timescale) to the observed GODAS monthly mean ocean temperature below a depth of 10 m. Nudging was not applied in the upper 10 m. The time step of SIT and AGCMs exchange ocean surface fluxes varying associated with the model resolution, which is 720, 1800, and 900 s in ECHAM-SIT, CAM5-SIT, and HiRAM-SIT, respectively. AGCMs were coupled with the SIT in the tropical region between 30^∘ S and 30^∘ N and forced by prescribed monthly mean OISST outside this tropical belt.

The experiments comprised three sets of coupled AGCM simulations (ECHAM5-SIT, CAM5-SIT, and HiRAM-SIT) and standalone AGCM simulations forced by observed monthly mean OISST (ECHAM5, CAM5, and HiRAM) from 1985 to 2005. The experiments were designed to evaluate the effect of atmosphere–ocean coupling on MJO simulations. Table 1 presents the model and experiment details.

Table 1 Detailed information about the models and experiments.

2.3 Methodology

The analysis focused on the boreal cool season (November–April), which is when the eastward propagation tendency of the MJO is the most prominent. We used the CLIVAR MJO Working Group diagnostics package (CLIVAR, 2009) and a 20–100 d filter to analyze intraseasonal variability. The MJO phase composites were computed using the real-time multivariate MJO index (Wheeler and Hendon, 2004), defined as the leading pair of principal components of intraseasonal OLR, and 850 and 200 hPa zonal winds in the tropics.

The vertically integrated moist static energy (MSE) budget was diagnosed based on the following equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{〈{\displaystyle \frac{\partial h}{\partial t}}〉}^{\prime }& =-{〈u{\displaystyle \frac{\partial h}{\partial x}}〉}^{\prime }-{〈v{\displaystyle \frac{\partial h}{\partial y}}〉}^{\prime }-{〈\mathit{\omega }{\displaystyle \frac{\partial h}{\partial p}}〉}^{\prime }+{〈\mathrm{LW}〉}^{\prime }\\ & +{〈\mathrm{SW}〉}^{\prime }+{〈\mathrm{SH}〉}^{\prime }+\langle \mathrm{LH}{\rangle }^{\prime }\end{array},\end{array}$$

where h is the MSE ($h=cpT+gz+Lq$); u and v are the zonal and meridional velocities, respectively; ω is the vertical pressure velocity; LW and SW are the longwave and shortwave radiation fluxes, respectively; and LH and SH are the latent and sensible surface heat fluxes, respectively. The mass-weighted vertical integration from the surface to 200 hPa is denoted as 〈⋅〉, and the intraseasonal anomalies are represented as 〈⋅〉', which were isolated using a 20–100 d bandpass Lanczos filter (Duchon, 1979).

3 Results

3.1 MJO simulations: ECHAM5-SIT, CAM5-SIT, and HiRAM-SIT

3.1.1 General structure

We compared simulated MJO characteristics using three coupled and uncoupled AGCMs. Figure 1 shows the wavenumber–frequency spectra of simulated 850 hPa zonal wind (shading) and precipitation (contours). All three uncoupled AGCMs (hereafter referred to as AGCMs) simulated intraseasonal signals with a lower frequency than was observed and overestimated the westward propagation with periods greater than 80 d (Fig. 1e–g). ECHAM5 and HiRAM simulated signals of wavenumbers 1–3 instead of the observed wavenumber 1 in 850 hPa zonal wind. These results show that all three AGCMs simulated stationary fluctuations with low frequency that were not consistent with the observations. In contrast, coupled AGCMs realistically reproduce the observed spectral characteristics and strength of the eastward propagation at wavenumbers 1 to 2 in 850 hPa zonal wind (Fig. 1b–d). Although HiRAM simulated eastward propagation in a wider frequency spectrum than observed, the coupled model clearly displays improvements in the MJO simulation compared with the stationary intraseasonal fluctuation in the uncoupled simulation. Hovmöller diagrams presented in Fig. 2 illustrate the temporal evolution of 850 hPa zonal wind and precipitation in the tropics in both observations and simulations. All three models simulated either stationary (CAM5 and HiRAM) or weak eastward-propagating (ECHAM5) signals in AGCMs, but they more realistically simulated the eastward propagation of the MJO in the coupled models. However, the propagation in the ECHAM5-SIT is still slightly slower than was observed. The improvement obtained in coupled models suggests that active ocean–atmosphere interaction is crucial for successful MJO simulation.

Figure 1 Wavenumber–frequency spectra for equatorial 850 hPa zonal wind (shading; m² s⁻²) and precipitation (contours; mm² d⁻²) over 10^∘ S–10^∘ N from (a) observations and simulations using the (b–d) coupled and (e–g) uncoupled AGCMs.

Figure 2 The 10^∘ S–10^∘ N averaged lag–longitude diagrams of intraseasonal precipitation (shading) and 10 m zonal wind (contour) correlated against precipitation (for 10^∘ S–5^∘ N, 120–150^∘ E) from (a) observations and simulations using the (b–d) coupled and (e–g) uncoupled AGCM. The contour interval is 0.1.

3.1.2 Atmospheric and oceanic profiles

The composite MJO life cycle, featuring intraseasonal OLR and 10 m surface wind anomalies for boreal winter in eight phases (following Wheeler and Hendon , 2004), is displayed in Fig. 3. All three coupled models simulated realistic MJO with enhanced circulations and propagation tendency compared to the uncoupled models. The MJO in phase 4, when deep convection is the strongest over the Maritime Continent (MC), demonstrates the large-scale zonally overturning circulation coupling with the convection (Fig. 4). The positive heating region in the coupled experiment is significantly enlarged, deepened, and tilted westward with increasing height compared to those in the uncoupled experiment. Correspondingly, the convective circulation envelope of the MJO is thicker and longitudinally wider in coupled experiments. The strong convection is associated with much enhanced low-level moisture convergence (green contours). Furthermore, the area of positive rainfall anomaly in the coupled experiment becomes larger, and the sea level pressure anomaly is meridionally more confined, exhibiting the characteristics of intensified Kelvin-wave-like perturbations to the east of the deep convection. This enhancement of low-level moisture convergence is consistent with the frictional wave-conditional instability of the second kind of mechanism (Frictional wave CISK; Wang and Rui, 1990; Kang et al., 2013). The enhancement of the Kelvin wave can be observed in the symmetric wavenumber–frequency spectra (Fig. 5). The spectra between 0 and 0.35 d⁻¹ are presented to highlight the MJO and equatorial Kelvin waves. The coherence at wavenumbers of 2–4 for the 10–20 d period is simulated to be stronger in the three coupled models than in the uncoupled models.

Figure 3 Composite November–April 20–100 d OLR (W m⁻¹; shading) and 10 m surface wind anomalies (m s⁻¹; vectors) as a function of the MJO phase in (a) ECHAM5-SIT, (b) ECHAM5, (c) CAM5-SIT, (d) CAM5, (e) HiRAM-SIT, and (f) HiRAM. Vectors <0.6 m s⁻¹ are not shown. The reference vector is shown (in units of m s⁻¹) at the bottom right. The number of days used to generate the composite for each phase is shown to the right of each panel.

Figure 4 Structure of simulated MJO in phase 4. The longitude–height cross sections (averaged over 10^∘ S–EQ) of the MJO-scaled wind circulation (vector, u: m s⁻¹; omega: 10⁻² Pa s⁻¹), Q1 (shading, unit: K d⁻¹), and the horizontal moisture convergence (green contour, unit: 10⁻⁶ g kg⁻¹ s⁻¹) from (a) observations and simulations using the (b–d) coupled and (e–g) uncoupled AGCMs. The contour interval of the moisture convergence is 8 × 10⁻⁶ g kg⁻¹ s⁻¹; the solid line indicates positive values. Precipitation (shading, unit: mm d⁻¹) and sea level pressure (contour, unit: hPa) are also shown. The contour interval of sea level pressure is 30 hPa; the dashed line indicates negative values.

Figure 5 Symmetric wavenumber–frequency spectra of 850 hPa zonal wind averaged over 10^∘ N–10^∘ S using the (a, c, e) coupled and (b, d, f) uncoupled AGCMs (in m² s⁻²).

In addition to the atmospheric structure, the SST (Fig. S1) and vertical profile of ocean temperature examined are presented in Fig. 6. The observed SST variation in MJO variability is reproduced well in all three coupled models (Fig. S1). The warm SST leads the main MJO convection by approximately 5–10 d, followed by the cold SST approximately 5–10 d later (Flatau et al., 1997; DeMott et al., 2015; Tseng et al., 2015). Moreover, the observed amplitude fluctuation (approximately 0.5 to 1 ^∘C) is realistically simulated. The observed ocean temperature profiles, characterized by the warm layer, along the Equator from the Indian Ocean to the western Pacific are simulated well in the three coupled models (Fig. 6). Meanwhile, simulated temperature anomalies are larger in ECHAM5-SIT than in CAM5-SIT and HiRAM-SIT. Figure S2 shows the fluctuations of observed SST and simulated SST in three sets of coupled and uncoupled models. There is no fluctuation (as expected) in uncoupled simulations, whereas the simulated SST fluctuates with phases similar to those observed at different locations. The amplitudes in ECHAM5-SIT and CAM5-SIT are similar to the observed values, whereas those in HiRAM-SIT seem to be smaller in the western Pacific. The differences between models are likely due to the different atmospheric model configurations because they were otherwise coupled to the same 1-D ocean model. Since the atmosphere is the main driver to extract heat from the ocean, different responses of atmospheric models likely have different effects on SST. Pinpointing the cause of quantitative differences between models will require further detailed analysis. The consistent results in all three coupled models support the conclusion of Tseng et al. (2015) that resolving fine vertical resolution in the upper ocean improves the simulation of the warm layer and MJO propagation and variability. Our results further demonstrate that the effect of atmosphere–ocean coupling on the MJO could be independent of AGCMs with different configurations and atmospheric physical parameterizations and that coupling seems to be a more fundamental approach.

Figure 6 Vertical ocean temperature (^∘C) profiles with respect to MJO phases for intraseasonal anomalies (i.e., with 20–100 d filtering) in observations and simulations using coupled models. Observations are in line with data from TAO. Because of storage limitations, only 3 and 10 m water temperatures are presented in the HiRAM-SIT simulation.

3.1.3 Performance comparison

Model performance is summarized in Fig. 7. The scatter plot shows the power ratio of east–west-propagating waves (x axis) versus the pattern correlation between the simulated and observed precipitation anomaly in Hovmöller diagrams (Fig. 2; y axis). The east:west ratio was calculated by dividing the eastward-propagating power by the westward-propagating power of 850 hPa zonal wind summed over wavenumbers of 1–2 and a period of 30–80 d. Compared with the observations, coupled simulations (marked by circles) exhibit better simulation than uncoupled simulations (marked by asterisks). A comparison of combined explained variance using real-time multivariate MJO series 1 (RMM1) and 2 (RMM2) (Fig. 7b) based on Wheeler and Hendon (2004) shows marked increases after coupling. The comparison demonstrates that coupling is essential for realistic MJO simulations.

Figure 7 Scatter plots of various MJO indices based on observations and experiments (Table 1). (a) The x axis is the power ratio of east–west-propagating waves. The east–west ratio was calculated by dividing the sum of eastward-propagating power by its westward-propagating counterpart within wavenumbers 1–3 (1–2 for zonal wind) and a period of 30–80 d. The y axis is the pattern correlation of precipitation and eastward propagation, as shown in Fig. 2. (b) Sum of RMM1 and RMM2 variances based on Wheeler and Hendon (2004).

3.2 Mechanism discussion

We applied the MSE budget to diagnose the moisture budget associated with the MJO. Figure 8 shows a Hovmöller diagram of MSE tendency averaged from the area between 10^∘ S and the Equator and the overlying precipitation anomalies. MSE tendency derived from reanalysis fluctuates in quadrature with precipitation anomaly with positive (negative) MSE tendency, leading (lagging) major convection by approximately one to two phases (DeMott et al., 2015, 2016, 2019). Coupled models simulate stronger eastward propagation in the MSE tendency and precipitation anomalies and realistic phase lag between the two. Stronger MSE tendencies in coupled simulations are observed in ECHAM5 and HiRAM but are less clear in CAM5. Figure 8d, g, and j show the differences between coupled and uncoupled simulations. One notable feature is the positive (negative) MSE tendency preceding positive (negative) precipitation anomaly, which preconditions environments for eastward propagation of active (inactive) convection and associated circulation. Next, we diagnosed the relative contribution of each term in Eq. (1) to the MSE tendency with a focus on the MC, where the largest positive MSE tendency and precipitation anomaly were found.

Figure 8 Hovmöller diagrams averaged over the area between 10^∘ S and the Equator for MSE (shading; J kg⁻¹) and precipitation (contour; mm d⁻¹) composites following the RMM index from (a) observations and simulations using the (b, e, j) coupled and (c, f, k) uncoupled AGCMs and (d, i, l) their difference. The contour interval is precipitation anomalies.

3.2.1 Preconditioning phase

Following the peak MSE tendency over the MC (120–150^∘ E) during phase 2 (Fig. 8d, g, and j), values of each term contributing to the column-integrated MSE tendency in Eq. (1) preceding the deep convection over the MC area (10^∘ S–0^∘ N/S, 120–150^∘ E) are shown in Fig. 9. Vertical advection is the dominant term with the major compensation from longwave radiation during phase 2 when convection is still in the eastern Indian Ocean, as identified by Wang and Li (2020). Moreover, the LH term is consistent within all three models and contributes less negative MSE tendency in coupled models than AGCMs. The results show that the contribution comes from the LH term in this early phase stage. The LH effect was overlooked in Tseng et al. (2015) because of the weak MJO variability in coupled simulations. However, this negative LH bias becomes one of the key factors in enhancing the leading MSE tendency during the MJO preconditioning phases. This suggests that the surface latent flux bias in AGCMs can be corrected by involving the coupling process in the preconditioning phase. Generally, coupling improves the budget simulation. The positive contribution of vertical advection and negative contribution of LH in MSE tendency is closer to being realistic in the coupled simulations during the initial phase of the MJO.

Figure 9 Model-simulated column-integrated MSE budget terms (J kg⁻¹ s⁻¹) during phase 2 of the MJO. Black, red, and blue represent the data from the observations, Nordeng scheme simulations, and Tiedtke scheme simulations, respectively. The averaged domain is 10^∘ S–0^∘ N/S, 120–150^∘ E.

Figure 10 Phase 4 of the column-integrated MSE tendency (shading; J kg⁻¹ s⁻¹) and precipitation (contours; mm d⁻¹) based on (a) observations, (b) ECHAM5-SIT, (c) ECHAM5, (d) CAM5-SIT, (e) CAM5, (g) HiRAM-SIT, and (f) HiRAM. The nine-point local smoothing is applied here to the intraseasonal precipitation variance of HiRAM (contours only).

3.2.2 Phase of strongest convection over the MC

We compared the spatial distribution of MSE and precipitation in phase 4 when convection was the strongest in the MC (Fig. 10). In the observations, the main convection occurs in the MC from 90 to 150^∘ E. A positive MSE tendency with a maximum value near 10^∘ N and 10^∘ S is identified in the east of the MJO convection centered near the Equator. Meanwhile, a negative integrated MSE tendency is found in the west of the MJO convection, and the meridionally confined structure near the Equator exhibits the characteristics of the equatorial Kelvin wave embedded in the MJO. Clearly, coupled models outperform uncoupled models in reproducing these signals. To quantify the contribution of coupling to the improvement, we follow Jiang et al. (2018) to project all MSE terms to the observations (Fig. 11). The dominant contribution of horizontal advection to the MSE tendency in the observations (Fig. 11a) is simulated well in the coupled simulations (but not in the uncoupled simulations) of ECHAM5 and CAM5 (Fig. 11b and c). Although a similar dominant effect was observed in both simulation types in HiRAM, it is enhanced in the coupled simulation (Fig. 11d). The horizontal advection term is further decomposed into zonal and meridional components (Fig. 11e–h); both components have a positive contribution, but the meridional component has a larger amplitude. Furthermore, the uncoupled ECHAM5 and CAM5 models simulate unrealistic features: positive contribution from zonal advection but negative contribution from meridional advection. In contrast, coupled models simulate the dominance of meridional advection well. In HiRAM, the uncoupled model simulates almost equally positive contributions from both terms. However, the coupled model simulates a larger contribution from meridional advection. We further decompose the meridional advection to assess the relative contributions of an intraseasonal perturbation and the mean state. Consistent with the observations (Fig. 11i), the meridional advection by intraseasonal flow ($-v^{\prime }\frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial y}$) is the main factor in improving the simulations in the coupled models (Fig. 11j–l). Our results are consistent with those of Jiang et al. (2018). To evaluate the relative contribution of intraseasonal circulation and background moisture, we further diagnosed changes in $\mathrm{\Delta }(-v^{\prime }\frac{\partial \stackrel{\mathrm{‾}}{h}}{\partial y})$ at phase 4. Here the overbar shows that the time mean and prime represents intraseasonal anomaly. Changes in the MJO meridional advection term for coupled experiments relative to uncoupled can be written as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{cc}\mathrm{\Delta }\left(-v^{\prime }\frac{\partial \stackrel{\mathrm{‾}}{m}}{\partial y}\right)=-\mathrm{\Delta }{v}^{\prime }{\left(\frac{\partial \stackrel{\mathrm{‾}}{m}}{\partial y}\right)}_{\mathrm{uncoupled}}& -{\left(v^{\prime }\right)}_{\mathrm{uncoupled}}\mathrm{\Delta }\left(\frac{\partial \stackrel{\mathrm{‾}}{m}}{\partial y}\right)\\ \left(\mathrm{a}\right)& \left(\mathrm{b}\right)\\ -\mathrm{\Delta }{v}^{\prime }\mathrm{\Delta }\left(\frac{\partial \stackrel{\mathrm{‾}}{m}}{\partial y}\right)& \\ \left(\mathrm{c}\right),& \end{array}\end{array}$$

where Δ represents the coupled–uncoupled change. The terms (a)–(c) in Eq. (2) are presented as bar charts in Fig. 12. The change in the intraseasonal circulation in the meridional circulation is the dominant factor in coupled simulations relative to uncoupled experiments. The instantaneous SST horizontal distribution dominates this moisture budget change due to the atmosphere–ocean coupling effect. Therefore, the change of varying moisture induces the intraseasonal circulation change. The results confirm that the dominance of dynamic influence over thermodynamic response to atmosphere–ocean coupling is essential in improving MJO simulations.

Figure 11 (a–d) The relative role of each MSE component of phase 4 through the projection of the spatial pattern of each MSE budget over the MC (domain) onto the total MSE tendency pattern (Fig. 8a). (e–h) Decomposite of the total horizontal MSE advection based on zonal and meridional components. (i–l) Relative role of the meridional horizontal MSE advection based on the MJO circulation and the mean state of moisture.

Figure 12 Bar chart of the relative contribution of intraseasonal convergence and background moisture between the coupled and uncoupled changes in MJO phase 4.

3.3 Discussion: mean state and intraseasonal variance

We examined the simulated mean state, which has been suggested a key factor affecting MJO simulations (Inness et al., 2003; Watterson and Syktus, 2007; Kim et al., 2009, 2011, 2014; Jiang et al., 2018, 2020). The three models exhibited different tropical SST responses to coupling (Fig. S3e). Over the warm pool area, CAM-SIT and HiRAM-SIT underestimate the SST, whereas ECHAM5-SIT overestimates the SST. Note that warm SST bias in the eastern tropical Pacific was simulated in the three models due to the lack of oceanic circulation in the SIT. The simulated zonal wind in the three models (Fig. S3b–d) demonstrated different responses to coupling. Figure S2c show the 850 hPa zonal wind differences between coupled and uncoupled models (shading) and the total field in uncoupled models (contours). Figure S3f–h show the 10^∘ S–0^∘ N/S averaged 850 hPa zonal wind in the coupled and uncoupled models. In ECHAM5-SIT, the westerly wind is slightly enhanced in the eastern Indian Ocean but decreases in the western Indian Ocean and western Pacific Ocean. In CAM5-SIT, westerly wind reduces in the Indian Ocean but enhances over the western Pacific. HiRAM-SIT has similar changes to those in ECHAM5-SIT, i.e., it decreases over the MC area but increases in the western Indian Ocean and Pacific. Generally, the three models disagree on the zonal wind mean state changes in response to coupling.

The mean moisture changes are substantially enhanced over the tropical areas in ECHAM5 after coupling (Fig. S4b and e). However, in CAM5 and HiRAM no clear change was observed to the south of the Equator, but strong drying was observed to the north for the same models (Fig. S4c, d, f, and g). The only common feature among the three models is the enhanced meridional gradient of mean moisture, which is consistent with many previous studies (Kim et al., 2014; Jiang et al., 2018; Ahn et al., 2020). Our budget analysis demonstrated that the meridional transport by the intraseasonal meridional circulation is the dominant term. It also showed that the meridional gradient of mean moisture is the secondary effect in enhancing MJO simulations by coupling. After coupling, the mean precipitation changes are more consistent among the three models (Fig. S5). One of the major changes is the southward shift of the major precipitation zone, resulting in precipitation increases over the regions south of the Equator (except in the MC). Similarly, the precipitation intraseasonal variance (20–100 d filtered) was markedly enhanced in these regions (Fig. S6). ECHAM5-SIT exhibits a relatively minor increase over the western MC. In contrast, HiRAM-SIT exhibits the strongest enhancement, particularly in the Indian Ocean. Generally, all three coupled models enhance the intraseasonal signals over the tropics with discrepancies in their level of detail. Meanwhile, the model mean state does not substantially improve after coupling. Thus, in this study, the mean state is not the main contribution to enhancing the MJO simulation after coupling. Instead, coupling leading to rigorous atmosphere–ocean interaction on intraseasonal timescales is likely the reason for the improving MJO simulation.

4 Conclusion

This study used a one-column, TKE-type ocean mixed-layer model (SIT) coupled with AGCMs to improve MJO simulation. SIT is designed to have fine layers near the surface that simulate the warm layer, cool skin, and diurnal fluctuations well. This refined discretization under the ocean surface in SIT provides improved SST simulation, which provides a more realistic air–sea interaction. Coupling SIT with ECHAM5, CAM5, and HiRAM significantly improves the MJO simulation in the three AGCMs compared with that of the prescribed SST-driven AGCMs. The vertical cross section indicates that the strengthened low-level convergence during the preconditioning phase is better simulated in the coupled experiment. Furthermore, the phase variation and amplitude of the SST and ocean temperature under the surface can be realistically simulated. Our results reveal that the MJO can be realistically simulated in terms of strength, period, and propagation speed by increasing the vertical resolution of the one-column ocean model to better resolve the upper-ocean warm layer.

The MSE budget analysis revealed that the coupling effects during the preconditioning and mature phases exhibit different contributions. During the preconditioning phase, the positive contribution of vertical advection and the negative contribution of LH in MSE tendency are closer to realistic values in coupled simulations during the initial phase of the MJO. Additionally, the meridional component of the horizontal advection term is the dominant term during the mature phase of the strongest convection in the MC, enhancing the simulation after coupling. Improved meridional circulation is essential in the coupled simulations that outperformed uncoupled experiments. The results confirm that the dominance of dynamic influence over thermodynamic influence in response to the atmosphere–ocean coupling is the key process in improving MJO simulations. In summary, this study suggests two major enhancements of the coupling process. First, the underestimated surface LH bias in AGCMs can be corrected during the preconditioning phase of the MJO over the MC. Second, during the strongest convection phase over the MC, the change in intraseasonal circulation in the meridional circulation is the dominant factor in coupled simulations relative to uncoupled experiments. Although many studies have indicated the key role played by the mean state, the mean state in our simulations provides only a secondary contribution to enhancing the MJO simulation, with coupling being the main contributor. For example, zonal wind and precipitation changed inconsistently in the three models after coupling. Instead, the meridional gradient of the mean moisture and intraseasonal precipitation variance has a better relationship after coupling. Therefore, coupling leading to rigorous atmosphere–ocean interaction in the intraseasonal time scale, but no change in mean states, is likely the reason for MJO simulation improvement. This study supports previous findings (Tseng et al., 2015) that show that enhancing atmosphere–ocean coupling by considering an extremely high vertical resolution in the first few meters of the ocean model improves MJO simulations. It also supports that this improvement is independent of AGCMs with different configurations and physical parameterization schemes. Resolving the atmosphere–ocean coupling may be more beneficial than modifying the atmospheric physical parameterization schemes in general circulation models. In brief, this study suggested the effectiveness of air–sea coupling for improving MJO simulation in a climate model and demonstrated the critical effect that being able to simulate the warm layer has on the results. Additionally, the findings presented here enhance our understanding of the physical processes that shape the characteristics of the MJO.