2025 Volume 103 Issue 6 Pages 679-706
We investigated the role of air-sea interaction in the intra-seasonal fluctuation of the monsoon trough (MT) over the western North Pacific (WNP) during boreal summer by comparing a global non-hydrostatic atmospheric numerical model “NICAM” and its ocean coupled version “NICOCO.” We conducted a series of 10-member ensemble experiments from 10 July to 1 September 2020, with NICOCO, NICAM forced by sea surface temperature (SST) that only fluctuates along the climatological seasonal cycle (NICAM-A), and NICAM forced by dailymean SST simulated in the NICOCO experiments (NICAM-N). All models simulated small eastward extension of the MT, characteristic under a La Niña condition. NICAM-N slightly overestimated MT activities. Compared to NICOCO, NICAM-A significantly underestimated the northward propagation speed of the convective envelope associated with the MT. This difference was attributed to variations in SST fluctuations. Clouds associated with the MT reflect downward shortwave radiation, reducing ocean heating. Additionally, when convective activities intensify over the WNP, cyclonic circulations associated with convections accelerate (decelerate) the background westerly/southwesterly on the southern (northern) side of convections, enhancing (reducing) latent heat release. Consequently, NICOCO simulated dominant negative net downward surface heat flux and resultant negative SST anomaly on the southern side, whereas NICAM-A did not simulate such fluctuations since the effect of the atmosphere to the ocean is not considered. The meridionally asymmetric SST anomaly structure simulated in NICOCO likely supports the northward propagation of convections. These results suggest that air-sea interaction is important in predicting the fluctuation of the MT and convections over the WNP.
In Asia during boreal summer, south-westerly monsoonal wind dominates and the zonally elongated confluence region from the northern Indian Ocean to the western North Pacific (WNP), called the monsoon trough (MT), is generated. The surrounding area of the MT features a large meridional gradient of zonal winds and relative vorticity, leading to intense convective activity (Holland 1995). The MT over the WNP influences various weather phenomena over east/southeast Asia (e.g., Ritchie and Holland 1999; Yoshida and Ishikawa 2013; Huang et al. 2022). Ritchie and Holland (1999) and Yoshida and Ishikawa (2013) reported that 60–70 % of tropical cyclones over the WNP are affected by the meridional shear and the confluence of low-level zonal winds. Huang et al. (2022) also noted that the MT over the WNP caused extreme rainfall in Henan in July 2021. Therefore, deepening our understanding of the MT over the WNP is crucial for comprehending and predicting various phenomena over east and southeast Asia.
Climatologically, the MT onsets over the South China Sea in mid-May, extends eastward toward the Philippines Sea by late-July, and retreats in late October (Ueda and Yasunari 1996; Wu and Wang 2001; Wu 2002). Its activity differs depending on mean sea surface temperature (SST) over the tropical Pacific in a given year (Wu and Wang 2000; Wu et al. 2012), particularly in relation to the condition of El Niño-Southern Oscillation (ENSO). For instance, during El Niño conditions, when SST anomalies over the equatorial eastern Pacific are positive, the subtropical high shifts eastward (Wang et al. 2000), allowing the MT to extend eastward (Lau and Chan 1986; Wang and Wu 1997). The MT also exhibits intra-seasonal fluctuations (e.g., Ueda and Yasunari 1996; Fig. 1), which are related to certain types of intra-seasonal oscillations (ISOs) (Lau and Chan 1986; Zhou and Chan 2005). During boreal summer, boreal summer ISO (BSISO; Yasunari 1979) dominates in the tropics of the Northern Hemisphere. BSISO has a 30–60 day periodicity and propagates northward over the northern Indian Ocean and WNP, as well eastward along the equator.
However, these various variabilities of the MT have not been fully understood because the relevant phenomena are not only interacted nonlinearly but also have internal variabilities. Specifically, there is no consensus on the mechanism behind the northward propagation of BSISO, despite extensive discussion in previous studies (Kikuchi 2021). For example, Jiang et al. (2004) indicated that mean easterly vertical wind shear (VWS), characteristic over the Asian summer monsoon region during boreal summer, induces barotropic vorticity to the north of BSISO convection, which then propagates northward. The contribution of mean meridional VWS has also been suggested (Bellon and Sobel 2008; Dixit and Srinivasan 2011). Additionally, the roles of moisture advection in the planetary boundary layer (PBL) (Jiang et al. 2004) and air-sea interaction (Wang et al. 2018; Yang et al. 2020) have been highlighted.
Over the WNP during boreal summer, it is important to consider air-sea interaction when simulating convections (Wang et al. 2005; Wu and Kirtman 2007). Wu and Kirtman (2007) indicated that atmospheric forcing dominates in the tropical Indo-western Pacific region during boreal summer, and atmosphere-only numerical models perform poorly when simulating convections in this region. However, how SST, affected by the atmosphere, feedback to convections remains unclear. Indeed, the impact of air-sea interaction on the northward propagation of BSISO convections varies among numerical models (Neena et al. 2017). Therefore, the objective of this study is to understand how SST, influenced by MT fluctuations, feedback to the MT.
Many previous studies have examined the fluctuation of convections associated with the MT using observational and reanalysis data (e.g., Wu 2002; Wang et al. 2018). Wu (2002) investigated the relation between northeastward extension of the MT and the northeastward propagation of the SST warm region, and indicated that the MT-induced SST variability causes the eastward extension of the MT through cloud-radiation and convection-evaporation feedbacks. Wang et al. (2018) examined the role of air-sea interaction in the northward propagation of BSISO over the WNP, and quantitatively revealed that the preceding positive SST anomaly on the northern side of BSISO convection helps the convection propagate northward. These observational studies suggest that SST fluctuations influenced by the atmosphere play a significant role in the intra-seasonal fluctuations of convections over the WNP.
There are also previous studies that examined the role of air-sea interaction in the fluctuation of convections associated with the MT using numerical models (e.g., Fu et al. 2002; Fu and Wang 2004; Neena et al. 2017). Fu et al. (2002) analyzed the local and remote air-sea interactions in the tropical Indian and Pacific Oceans through sensitivity experiments, and suggested that air-sea coupling over the Indian Ocean and WNP modifies the representation of the Asian summer monsoon climatology. Fu and Wang (2004) compared an atmospheric-only general circulation model (AGCM) and ocean-coupled GCM (CGCM). They demonstrated that CGCM can simulate stronger BSISO than AGCM because SST fluctuations influenced by the atmosphere are considered in CGCM. However, they also suggested that air-sea interaction does not play a significant role in the period and phase speed of northward propagation of BSISO. Neena et al. (2017) evaluated the representation of BSISO in 20-year climate simulations with 22 AGCMs and 5 CGCMs. This study concluded that air-sea coupling is not necessary for BSISO convections to propagate northward. These results highlight the uncertainty regarding the role of air-sea interaction in simulating convections over the WNP, indicating that there are still issues to be addressed.
Most previous studies investigating convections associated with the Asian summer monsoon using GCMs relied on cumulus and/or deep convective parameterizations due to coarse horizontal resolution (e.g., Fu et al. 2002; Fu and Wang 2003; Neena et al. 2017; Yang et al. 2020). These parameterizations often introduce large uncertainties, impacting the simulation of convection which characterizes the MT. Additionally, such GCMs have biases including the unrealistic location of the VWS in the background field (e.g., Fu and Wang 2004). Given that the MT affects the geneses of tropical cyclones, which can cause SST cooling through ocean dynamical effect (i.e., upwelling) (Price 1981), high horizontal resolution is required to simulate tropical cyclones. To address these challenges, we used a global non-hydrostatic numerical model “NICAM” (Tomita and Satoh 2004; Satoh et al. 2008, 2014), and its ocean coupled version “NICOCO” (Miyakawa et al. 2017). Due to its high horizonal resolution, NICAM is capable of representing convection without relying on convection parameterization. NICAM is also outstanding in simulating tropical phenomena such as BSISO and tropical cyclones (e.g., Shibuya et al. 2021; Yamada et al. 2023), particularly with regard to clouds and precipitation. Furthermore, recent improvement has reduced model biases, including the double Intertropical Convergence Zone, allowing NICAM to realistically simulate the Asian summer monsoon (Takasuka et al. 2024). Thus, NICAM and NICOCO have advantages for examining the MT fluctuations in terms of air-sea interaction. To the best of our knowledge, there are no studies that have examined the role of air-sea interaction in the intra-seasonal fluctuation of the MT with not only high-resolved but also air-sea coupled models.
This study examines the role of air-sea coupling in the intra-seasonal fluctuation of the MT over the WNP. Section 2 describes the numerical models, data, and experiment settings. Model performance and differences among experiments are presented in Section 3. Section 4 discusses our findings, and the last section is devoted to conclusions.
In this study, we used the Non-hydrostatic Icosahedral Atmospheric Model (NICAM; Tomita and Satoh 2004; Satoh et al. 2008, 2014) and its ocean-coupled version “NICOCO” (Miyakawa et al. 2017), which is combined with the CCSR Ocean Component Model (COCO; Hasumi 2006) via “J-Cup” (Arakawa et al., 2011, 2020). NICAM and COCO were coupled every 30 minutes. NICAM was configured with an approximately 14 km horizontal grid spacing and 78 vertical layers (with a model top at about 50 km). COCO was configured with about 0.25 degrees horizontal grid spacing and 63 vertical levels (with the deepest model-bottom at about 7300 m). The ocean was initialized with a 10-day spin-up from reanalysis data to represent a realistic ocean field, as described in Section 2.4. Additionally, a simple flux correction (Masunaga et al. 2023) was applied in the NICOCO experiments to suppress SST drift caused by heat imbalance between NICAM and COCO.
The physics schemes used in NICAM are the same as those in Takasuka et al. (2024). The cloud microphysics scheme is an updated version of NICAM Single-Moment Water 6 (NSW6; Tomita 2008), which is based on satellite observations (Roh and Satoh 2014; Roh et al. 2017). Turbulence is calculated using the Mellor-Yamada-Nakanishi-Niino level2 (MYNN2; Nakanishi and Niino 2006; Noda et al. 2010) and the Leonard term (Leonard 1975; Germano 1986). Radiation is handled by the Model Simulation radiation TRaNsfer code X (MSTRN-X; Sekiguchi and Nakajima 2008). The land surface is modeled using the Minimal Advanced Treatments of Surface Interaction and RunOff (MATSIRO; Takata et al. 2003). Sea surface heat and momentum fluxes are computed using the Louis-type bulk scheme (Louis 1979; Fairall et al. 2003; Moon et al. 2007). Parameters related to cloud microphysics and the Leonard term follow the “MP-LEO1” configuration in Takasuka et al. (2024), which reasonably reproduces the MT. The physics schemes used in COCO are the same as those in Masunaga et al. (2023). The turbulence closure scheme (Noh and Kim 1999) is applied to the surface mixed layer. Second-order momentum conservation scheme (Prather 1986) is used for tracer advection. Bi-harmonic Smagorinsky-like viscosity (Griffies and Hallberg 2000) is used for horizontal viscosity. The model configurations are summarized in Tables 1 and 2.


For global ensemble experiments, we used the AFES-LETKF Experimental Ensemble Reanalysis (ALERA; Miyoshi et al. 2007) version3 (ALERA3; Yamazaki et al. 2023) for atmospheric initial conditions. ALERA is an atmospheric ensemble dataset generated through the Local Ensemble Transform Kalman Filter (LETKF; Hunt et al. 2007) data assimilation system for the Atmospheric General Circulation Model for the Earth Simulator (AFES). In ALERA3, observational data from the National Centers for Environment Prediction (NCEP) and brightness temperature data from the Advanced Microwave Sounding Unit-A (AMSU-A) are assimilated (Terasaki and Miyoshi 2017). The horizontal and vertical resolutions of ALERA3 are 1.25 degrees and 37 layers, respectively.
To initialize the ocean in COCO (Section 2.4), we used daily-mean ocean reanalysis data from the Global Ocean Physics Reanalysis (GLORYS12v1; Lellouche et al. 2018; Drévillon et al. 2022), with horizontal and vertical resolutions of 1/12 degrees and 50 levels (with the deepest-bottom at about 5700 m). This dataset assimilates data from the ocean model NEMO (Madec and Imbard 1996) and observational data using a reduced-order Kalman filter (Brasseur and Verron 2006). The upper boundary condition is provided by ERA-Interim (Dee et al. 2011) or ERA5 (Hersbach et al. 2020). Observational data include sea surface height, SST, and ice concentration from satellite observations, as well as vertical profiles of temperature and salinity from in-situ observations. Although there are biases in the slow evolution of temperature and salinity, these biases are corrected using a 3-dimensional variational scheme (Lellouche et al. 2013).
For COCO experiments requiring ocean initialization with reanalysis data (details in Section 2.4), the upper boundary condition is provided by the Japanese 55-year atmospheric reanalysis designed for driving ocean-sea ice models (JRA55-do; Tsujino et al. 2018). JRA55-do, created by adjusting the JRA55 (Kobayashi et al. 2015) atmospheric reanalysis with satellite data and other atmospheric reanalysis datasets, is used specifically for ocean modeling.
ERA5 data are also used to create the bottom boundary condition in NICAM experiments and evaluate the realism of the experiments. Additionally, outgoing longwave radiation (OLR) data (Liebmann and Smith 1996) and SST data (Huang et al. 2021), which are provided by National Oceanic and Atmospheric Administration (NOAA), are used for the same purpose.
2.3 Experiment settingsWe conducted 10-member ensemble numerical experiments with NICAM and NICOCO from 10 July to 1 September 2020. In NICAM experiments, two types of SST were prescribed for the bottom boundary condition: one was the daily-mean SST that fluctuates along the climatological seasonal cycle (hereinafter SST-A), and the other was the daily-mean SST simulated by NICOCO (hereinafter SST-N). To compare the effect of biases, ERA5 and NOAA data were used for SST-A. In this study, we refer to the NICAM experiments forced by SST-A with ERA5 data as “NICAM-A (ERA5),” those with NOAA data as “NICAM-A (NOAA),” and those forced by SST-N as “NICAM-N.” Details of SST-A and SST-N are described below, and the experiment settings are summarized in Table 3.

SST-A was calculated as the sum of the climatology and an anomaly, both of which were derived from ERA5 and NOAA data. The climatology was calculated as 40-year mean over the period from 1980 to 2019. The anomaly was defined, following the procedure in Miyakawa et al. (2014), as the difference between the climatology and the average from 3 to 9 July 2020, which reflects the La Niña condition in the 2020 summer. In a NICAM-N experiment, the SST-N distribution was obtained from a NICOCO experiment that used the same ALERA3 realization (i.e., ensemble number) with the NICAM-N experiment for the atmospheric initial condition. Specifically, the SST simulated in NICOCO using member-1 of ALERA3 was used as the bottom boundary condition in NICAM using member-1 of ALERA3, and so on for the other members. The primary difference between NICOCO and NICAM-N is the representation of diurnal variations depending on coupling and uncoupling. Such variations may have effects in simulating convections, which would be successfully simulated in a high-resolution numerical model.
2.4 Description of initializationBefore describing how to initialize the ocean in COCO, we explain the commonly used method for providing ocean initial data. COCO is the ocean component of a climate model, the Model for Interdisciplinary Research on Climate (MIROC; Tatebe et al. 2019), which focuses on the long-term climate variability and prediction of global warming. In climate experiments, the ocean initial data are generated by running COCO for a few hundred years with JRA55-do (referred to as “Method A” and schematically shown in Fig. 1a) to simulate deep ocean circulation. However, the ocean climatology for the model will have bias, and the state of the internal variabilities may also differ from reality when NICAM and COCO are coupled using Method A. Because standard deviation of SST around the Philippines, where the MT is located, exceeds about 0.1 °C (Wang et al. 2018; Fig. 6), it is desirable that the deviation from reality is less than 0.1 °C when conducting forecast experiment. Providing realistic initial ocean internal structure near the sea surface, such as mixed layer, is also crucial because the shallow layers are related to the ocean response to the atmosphere forcing. As shown in Section 3 (Figs. 3c, d), the ocean initial state differs largely from reality. To start experiments with ocean state closer to reality, we initialized the ocean with the ocean reanalysis data GLORYS12v1 (referred to as “Method B” and schematically shown in Fig. 1b).

Schematics describing the methods for creating the ocean initial field in NICOCO. (a) Method A; the ocean field after long-term simulation is used as the initial condition. (b) Method B; the initial condition in COCO is initialized with ocean reanalysis data, and after a spin-up period of approximately 10 days, the ocean field is provided as the initial condition.

Time series of sea surface height tendency. The black, blue, and red thick lines represent the global, tropical 30°S–30°N, and WNP region [100–180°E, 0°–30°N], respectively, using Method B. The black dotted line represents the global region using Method A.
Method B involves two steps. First, we interpolated temperature and salinity from GLORYS12v1 onto the COCO grid. Flow velocity, flux, and sea surface height were set to zero to avoid numerical instability caused by the geometric differences between the GLORYS12v1 and COCO grids. Sea surface height in COCO is defined as the bottom of sea ice in order to represent the dynamical flatness taking isostasy into consideration. Because specific weight of sea ice (snow) relative to seawater is 0.9 (0.3), sea surface height is calculated as
![]() |
where TI is the thickness of sea ice, TS is the thickness of snow over sea ice, AI is ice concentration, and the subscript i in each variable means the category based on sea ice thickness. Variables related to sea ice are classified into five categories according to sea ice thickness, with the minimum thickness for each category being 30, 60, 100, 250, and 500 cm, respectively. For ice-related variables, we used data from Method A since no corresponding data for ice internal temperature was available. The difference in sea ice between GLORYS12v1 and COCO is unlikely to affect the tropical or WNP region focused on in this study at a sub-seasonal timescale. Second, since no initial flow is present in the ocean field, COCO was run until geostrophic flows developed. This process is referred to as “geostrophic spin-up” in this study. During the geostrophic spin-up, the ocean was forced with temperature and salinity from the initialization time to prevent the ocean in COCO from diverging from the ocean reanalysis data. The restoring time scale for the body forcing was set to 72 hours, and the upper boundary condition was provided by JRA55-do, similar to Method A. The spun-up ocean field is then provided to NICOCO as the ocean initial data.
We calculated the tendency of sea surface height to assess the development of large-scale flows. Figure 2 shows the tendency in the global (black line), tropical (blue line), and WNP regions (red line). The tropical region is defined as [0°–360°E, 20°S–20°N] and the WNP region as [100–180°E, 0°–30°N]. The black dotted line in Fig. 2 shows the global sea surface height tendency using Method A, which varies no more than 1.0 cm h−1. In contrast, the black solid line shows that the tendency in Method B decreases rapidly and reaches values similar to Method A after approximately 8 days. This pattern is also shown in the tropical and WNP regions. These results suggest that, even if flow velocity and sea surface height are initially set to zero, providing temperature and salinity allows sea surface height to become quasi-stationary and geostrophic flow to develop within about 8 days.

Temperature at the initial time of numerical experiments. (a) and (b) show SST and temperature at 100 m depth from GLORYS12v1, respectively. (c) and (d) show the difference between COCO with Method A and GLORYS12v1. (e) and (f) show the difference between COCO with Method B and GLORYS12v1. Gray (White) regions indicate land regions at surface (100 m depth).
In this subsection, we check the reproducibility of the ocean states in NICOCO, and grasp whether the initialization with ocean reanalysis (Section 2.4) is controlling the ocean behavior as intended. Figure 3 illustrates SST and temperature at 100 m depth on 10 July 2020, the initial time for the NICOCO experiments. Figure 3a highlights that SST around Philippines is higher compared to other regions. This warm region is accurately represented in COCO using Method B (Fig. 3e), whereas it is underestimated in COCO with Method A (Fig. 3c). Similarly, Fig. 3f demonstrates that temperature at 100 m depth aligns more closely with observational data when initialized with Method B, in contrast to Method A (Fig. 3d). Salinity also shows improved realism with Method B (not shown). Therefore, the thermohaline structure at the initial time is better represented with Method B than with Method A.

Distributions of the flow velocity over the WNP on 10 July 2020. Left panels (a c, e) are at sea surface and right panels (b d, f) are at 100 m depth. (a, b) show GLORYS12v1, (c, d) show the results with Method A, and (e, f) show results with Method B. Gray (white) regions indicate land regions at surface (100 m depth).

Zonal-depth cross sections of temperature at (a, b) 5°N, (c, d) 15°N, and (e, f) 25°N. Left panels are GLORYS12v1 and right panels are the results from NICOCO member-1. Land regions are masked with white.

50-day mean SST. (a) is GLORYS12v1 data, (b) is the result from NICOCO, (c) is that from NICAM-A (ERA5), (d) is that from NICAM-A (NOAA). Gray regions indicate land regions at surface.
Figure 4 presents the flow velocity at the surface and at 100 m depth over the WNP for GLORYS12v1 and COCO with Method A and B. The flow field simulated by COCO using Method B more closely resembles that of GLORYS12v1 than that using Method A. However, velocities near the equator, south of 5°N, are underestimated, where the flow is not dominated by geostrophic balance. From Figs. 3 and 4, initializing the ocean with Method B enables more realistic numerical experiments with an improved oceanic initial field. It should be noted that circulations at deeper depths are not fully established, but their impacts on the intra-seasonal phenomena studied here are expected to be negligible.
In this study, we applied the flux adjustment method (Masunaga et al. 2023) to mitigate SST drift. The flux adjustment may affect the internal thermohaline structure, particularly the ocean mixed layer. Hence, we examined vertical structure of temperature. Figure 5 shows zonal-depth cross sections of temperature from GLORYS12v1 and one member (member-1) at 5, 15, and 25°N. These figures indicate that differences of temperature structures between GLORYS12v1 and model outputs are small, confirming that NICOCO effectively simulates the ocean mixed layer. However, the bias of the mixed layer depth is observed in the east of approximately 160°E near the equator (Fig. 5b), possibly due to underestimated flow velocities (Fig. 4). It should be noted that other members exhibit similar structures (not shown). To summarize this section, the ocean internal structure around Philippines throughout the calculation period was successfully simulated by initialization with ocean reanalysis data and flux adjustment method.
3.2 Background fieldsThis subsection focuses on the overall reproducibility of 50-day mean fields. Figure 6 illustrates 50-day mean SST. Observational data (Fig. 6a) shows that SST around the Philippines is high, exceeding 30 °C. NICOCO (Fig. 6b) captures the warm SST region more realistically than NICAM-A (ERA5) (Fig. 6c). It should be noted that total downward heat flux in NICAM is larger than in COCO over the WNP in summer 2020, and SST experiences a warming drift if the flux adjustment (Masunaga et al. 2023) is not applied (not shown). NICAM-A (NOAA) also simulates the warm region although underestimation is observed in the east of Philippines. SST-A (ERA5) underestimates the warm SST region by approximately 0.5 °C because SST from ERA5 is lower than that from NOAA over the WNP during this period.
Figure 7 shows mean OLR and horizontal winds at 850 hPa over the calculation period. In summer 2020, the MT extends less eastward and the confluence region of low-level zonal winds is located around Philippines, which is characteristic under La Niña condition (Fig. 7a). All experiments successfully simulate the structure although NICAM-N exhibits slightly stronger westerly (Figs. 7b–e). Regarding the MT, it is observed that OLR is low around Philippines (Fig 7a). However, each numerical experiment substantially overestimates OLR around Philippines (Figs. 7b–e), indicating that convective activity is suppressed. In particular, OLR in NICAM-A (ERA5) (Fig. 7c) is overestimated likely because SST is underestimated (Fig. 6c). We have also calculated the seasonal evolution of the Asian summer monsoon (Fig. S1) following Takahashi et al. (2020), and confirmed that the monsoon in ERA5 and each numerical experiment is in the mature stage during the simulation period.

50-day mean (shaded) OLR and (vector) horizontal wind at 850 hPa for (a) NOAA and ERA5, (b) NICOCO, (c) NICAM-A (ERA5), (d) NICAM-A (NOAA), and (e) NICAM-N. (b–e) are the results of all member ensemble mean.
Figure 8 shows mean VWS between 200 and 850 hPa. ERA5 exhibits dominant zonal VWS (Fig. 8a), which is common structure over the Asian monsoon region in boreal summer (e.g., Kikuchi 2021). As described in Section 1, some previous studies emphasized that mean zonal VWS contributes to the northward propagation of BSISO (Wang and Xie 1997; Jiang et al. 2004), which influences fluctuations of the MT (Zhou and Chan 2005). Because similar VWS structure is also simulated in all numerical experiments (Figs. 8b–e), VWS are unlikely to be dominant factor causing the differences described in the following subsections.

50-day mean VWS between 200 hPa and 850 hPa for (a) ERA5, (b) NICOCO, (c) NICAM-A (ERA5), (d) NICAM-A (NOAA), and (e) NICAM-N. (b–e) are the results of all member ensemble mean. Shading represents the magnitude of the shear, and white regions in (b–e) indicate terrain at 850 hPa.
As described in the previous subsection, the MT extends zonally, on average, to the Philippines during the calculation period (Fig. 7). However, the field is not steady, and westerlies and convections fluctuate. Figure 9 shows snapshots of OLR and horizontal winds at 850 hPa in NOAA and ERA5 data. On 5 August (Fig. 9a), the westerly extends to the east of the Philippines (about 130°E), and the confluence region and convection are observed around that region. The westerly and convection moves northward after 3 days (Fig. 9b), and then the extension of the westerly and convection around the Philippines are weakened on 14 August (Fig. 9d). On 17 August (Fig. 9e), the westerly extends to the east of the Philippines again followed by generation of another convection, which also move northward after 3 days (Fig. 9f). These suggest that the MT fluctuates with variations of the westerly and convection within about 2 weeks, which is also observed in each numerical experiment (not shown). In the following, we examine such fluctuations before and after the intensification of the MT.

Snapshots of (shaded) OLR and (vectors) horizontal winds at 850 hPa in NOAA and ERA5 data.
The metrics for objectively evaluating fluctuations of the MT varies across studies (e.g., Wang et al. 2008; Feng et al. 2020), such as low-level zonal winds, low-level relative vorticity, and OLR. In this study, we used the maximum of the westerly at 850 hPa around the Philippines [110–130°E, 5–20°N] as the criterion of the time when the MT becomes the most active. After determining that time for observation and each numerical experiment (hereinafter lag = 0), composite analyses based on the time were conducted to examine the MT fluctuations. Figure 10 shows composites of zonal winds (shaded) and OLR (contour) at lag = 0, with hatched regions in Figs. 10c–e indicating that the difference from NICOCO is statistically significant at the 95 % confidence level based on Welch’s t-test. Figures 10a–e represent that the MT extends to the east of the Philippines in each composite. It is also found that the MTs in NICAM experiments extend more eastward than the NICOCO experiment.

Composites of zonal winds (shaded) and OLR (contour) at lag = 0 for (a) NOAA and ERA5 data, (b) NICOCO, (c) NICAM-A (ERA5), (d) NICAM-A (NOAA), and (e) NICAM-N. The composite is based on the time the westerly at 850 hPa around the Philippines [110–130°E, 5–20°N] is the strongest. Contours are drawn every 20 W m−2 less than or equal to 200 W m−2. Hatched regions in (c–e) indicate that the difference from NICOCO is statistically significant at the 95 % confidence level based on the Welch’s t-test. White regions in (b–e) indicate terrain at 850 hPa.
Figure 11 shows composites of SST (shaded), OLR (contour), and horizontal wind at 850 hPa (vector) anomalies, defined as the difference from 50-day mean of each member, for 4 days before and after the time when the MT is the most active (i.e., lag = −4 – +4 days) simulated in NICOCO. The monsoonal westerly around the Philippines (i.e., the MT over the WNP) starts to extend eastward at lag = −4, the westerly and convections become dominant at lag = 0, and then they move northward and weaken at lag = 4. In the next paragraph, we examine how the westerly and convections following the MT fluctuates in each experiment. Please note that details of SST are described in Section 3.4.

Composites of (shaded) SST, (contour) OLR, and (vectors) horizontal wind at 850 hPa anomalies simulated in NICOCO. Contour is drawn every −20 W m−2 below 0 W m−2. Gray regions indicate land regions at surface.
Figure 12 shows the same composites as Fig. 10 except that meridional mean over 5–20°N and time series for lag = −10 – +10 days are shown. It is noted that the strong westerly at lag = −10 – −5 days in observational data is related to a tropical cyclone. From Fig. 12, the eastward extension of the westerly is slightly stronger in NICAM experiments than NICOCO experiments, consistent with Fig. 10. In addition, the strong westerly region spread a few days longer in NICAM than in NICOCO. For instance, the period when the westerly exceeds 3 m s−1 on 120°E is 5 days in NICOCO, 7 days in NICAM-A (ERA5), 12 days in NICAM-A (NOAA), and 10 days in NICAM-N. From these results, it is suggested that the MT extends farther eastward and tends to be persistent in air-sea uncoupled simulations compared to air-sea coupled ones.

The same composites as Fig. 10 but for meridional mean over 5–20°N mean and for lag = −10 – +10 days. The horizontal and vertical axes represent longitude and lag-time (positive lag indicates leading time), respectively. Contours are drawn every 10 W m−2 less than or equal to 210 W m−2. Hatched regions in (c–e) indicate that the difference from NICOCO is statistically significant at the 95 % confidence level based on the Welch’s t-test.
Previous studies suggest that the MT over the WNP is influenced by intra-seasonal oscillations such as BSISO (Zhou and Chan 2005). Because convections associated with these oscillations tend to propagate northward over the WNP, it is hypothesized that differences of meridional fluctuations in the MT may also be observed. Figure 13 shows the same composites as Fig. 12 but for zonal mean over 110–130°E. The westerly region and convections initially occur in the south and propagate northward as the MT extends eastward in observational data (Fig. 13a). This northward propagation is also seen in numerical experiments (Figs. 13b–e). Focusing on Figs. 13b–d, the westerly region and convections following the MT in NICAM-A (ERA5) and NICAM-A (NOAA) propagate more slowly than in NICOCO. On the other hand, the westerly and convections following the MT in NICAM-N propagate almost as fast as those in NICOCO (Figs. 13b, e), although the westerly in NICAM-N is stronger.

The same composites as Fig. 12 but for zonal mean over 110–130°E. The horizontal and vertical axes indicate latitude and lag-time, respectively.
In the results above, the maximum of the westerly at 850 hPa around the Philippines was used as the composite metric. We also used other metrics and did the same composites to examine the sensitivity of the metric to results. The same composites as Fig. 13 when using the maximum of the westerly at 850 hPa over a narrower region [115–125°E, 10–15°N], for the sensitivity of the selected region, and the minimum of OLR around the Philippines, for the sensitivity of the selected variable, as the metric are shown in Fig. S2 and Fig. S3, respectively. In both cases, the slower fluctuation of the westerly in NICAM-A experiments than in NICOCO was calculated, which is similar to Fig. 13.
3.4 SST fluctuationAs described in the previous subsection, the zonal and meridional fluctuations of the MT differed among experimental settings. Particularly noteworthy was the result in NICAM-A, which simulated a slower or underestimated meridional fluctuation of the MT possibly because SST in NICAM-A fluctuates only climatologically. This section explores the discrepancies in SST and convection fluctuations among the experiments. Figure 14 presents lag composites of 110–130°E zonal mean of (a) SST, (b) SST tendency, (c) shortwave radiation, (d) longwave radiation, (e) latent heat flux, (f) sensible heat flux, and OLR (contour) anomalies, where anomaly is defined as the difference from 50-day mean of each member. Note that downward is positive for shortwave radiation, while upward is positive for longwave radiation and latent and sensible heat fluxes. It is noted that the 95 % confidence range based on a t-distribution for each variable is shown in Fig. S-4. Figures 14a and 14b indicate that SST decreases significantly around negative OLR anomaly regions (i.e., active convections) compared to relatively inactive convections. Around active convections, shortwave radiation anomaly is negative and corresponds exactly to negative OLR anomaly, as clouds following the development of convections reflect shortwave radiation that would otherwise reach the surface (Fig. 14c). On the other hand, latent heat flux anomaly is positive and larger on the southern side than the northern side of negative OLR anomaly regions (Fig. 14e). From Figs. 14d and 14f, longwave radiation and sensible heat flux anomalies are small, so they are unlikely not to contribute much to SST fluctuations.

Lag composites of (a) SST, (b) SST tendency, (c) surface downward shortwave heat flux, (d) surface upward longwave heat flux, (e) surface upward latent heat flux, (f) surface upward sensible heat flux, and (contour) OLR anomalies simulated in NICOCO. The composite is based on the same metric as Fig. 10. The horizontal and vertical axes represent latitude and lag, respectively. Contours in each figure are drawn every −20 W m−2 below 0 W m−2.
Figure 15 shows the composite analysis based on convection centers from lag = −5 day to lag = +5 day for ERA5 and each experiment. Convection centers are defined as the minimum OLR anomaly around the Philippines [110–130°E, 5–20°N], and the 110–130°E zonal mean is composited based on the latitude of the minimum location. Only values over the ocean are used for the composite, except for the OLR anomaly used for the composite, because we focus on the effect of surface heat flux fluctuations on SST. Longwave radiation and sensible heat flux are omitted in Fig. 15 since contributions of them to SST fluctuations are much smaller than shortwave radiation and latent heat flux (Fig. 14). Net downward surface heat flux averaged over every 5° latitude in Fig. 15 is summarized in Table 4, with the 95 % confidence ranges based on a t-distribution. From Fig. 15 and Table 4, net downward surface heat flux over the ocean (solid black line) is more negative on the southern portion than on the northern portion of convection. For instance, the value in the NICOCO experiment is −40.5 W m−2 over the −10° – −5° relative latitude range, while −10.2 W m−2 over the +5° – +10° range. The difference appears to increase as the distance from the center increases.

Composites based on the convection centers from lag = −5 day to lag = +5 day for (a) ERA5, (b) NICOCO, (c) NICAM-A (ERA5), (d) NICAM-A (NOAA), and (e) NICAM-N. The horizontal axis represents a relative latitude from convection centers, with the right (left) side of each figure indicating the northern (southern) side of the convection. In each figure, (shaded) SST, (black thick line) net downward surface heat flux, (black dashed line) OLR, (cyan line) surface shortwave radiation, (magenta line) surface latent heat flux, (orange line) SST tendency, and (green line) surface velocity. Whiskers drawn every 5 W m−2 indicate the 95 % confidence interval for each variable based on the t-distribution.

The difference can be explained as follows: OLR (black dashed line) and shortwave radiation flux (cyan line) values are nearly symmetric with respect to the convection because clouds reflect downward shortwave radiation. In contrast, latent heat flux (magenta line) is larger on the southern portion than on the northern portion because surface wind velocity (green line) is also meridional asymmetric, which is consistent with the fact that surface latent heat flux increases as surface wind velocities become large (Louis 1979). Figure 11 elucidates the meridional asymmetry of surface velocities. As convections around the Philippines become active, strong cyclonic circulations develop. In this region, where westerly/southwesterly winds following the MT generally blow (Fig. 7b), westerly/southwesterly winds on the southern portion of convections act to accelerate the background wind, while easterly/northeasterly winds on the northern portion of convections act to decelerate. Hence, cyclonic circulations generated around convections cause the meridional asymmetry of surface winds.
In the NICOCO experiments, SST tendency (orange line) is negative around convection centers, leading to a dominant negative SST anomaly region (blue shaded region) on the southern portion (Fig. 15b). However, the peak of SST tendency is located at +2° relative latitude, whereas the peak of net downward surface heat flux (solid black line) is centered on the convection. This northward shift in the SST tendency peak is attributed to the ocean upwelling induced by tropical cyclones (Price 1981) which is observed in four members. In fact, the confidence range of SST tendency at 15–25°N around lag = 0 day (Fig. S4) is large, which means variations of SST tendency in that region among ensemble members are large. On the other hand, NICAM experiments exhibit little SST fluctuations (orange lines in Figs. 15c–e) since the SST fluctuations simulated in NICOCO are not considered. Additionally, SST-A fluctuates only climatologically while SST-N is obtained from NICOCO simulations that include the atmospheric forcing simulated in NICOCO. Therefore, the difference between NICAM-A and NICAM-N is seen in the magnitude of SST anomaly. It is also noteworthy that net downward surface heat flux around convections is more negative in NICAM-A (ERA5) (average: −43.6 W m−2), NICAM-A (NOAA) (average: −51.4 W m−2), and NICAM-N (average: −48.7 W m−2) than in NICOCO (average: −41.2 W m−2) because more latent heat flux is released in NICAM due to larger surface winds (Fig. 15).
As shown in Section 3, the meridional fluctuation of the MT differs among experimental settings, particularly between NICOCO and NICAM-A (Fig. 13). This difference likely results from differences in SST fluctuations. In the NICOCO experiments, a meridional asymmetric structure in the SST anomaly, positive on the northern and negative on the southern side, is created because surface winds are stronger and larger latent heat flux exists over the southern portion of the MT, under the mean Asian monsoonal wind (Figs. 11, 15b). According to previous studies (e.g., Roxy and Tanimoto 2012; Ren et al. 2013; Wang et al. 2018), the development of convection is favored (suppressed) over positive (negative) SST anomalies. Thus, the asymmetric SST anomaly structure should be important for realistic northward propagation of convections (Fig. 13b). In the NICAM-A experiments, such SST fluctuations are not found, even though surface winds and latent heat flux are larger on the southern side of convections (Figs. 15c, d) since the prescribed SST fluctuates only climatologically. As consistent with the lack of the asymmetric SST anomaly structure, the northward propagation of convections and the westerly region following the MT are slow or underestimated (Figs. 13c, d), suggesting that air-sea coupling improves the speed of northward propagation of the positive SST anomaly, thereby improving the convections over the WNP. The mechanism described in this paragraph is consistent with DeMott et al. (2013), who indicated that SST anomaly contributes to the northward propagation of convections over the Western Pacific. The role is also observed in active/break cycle of monsoonal wind over the Bey of Bengal (Vecchi and Hendon 2002).
Such a role is similar to the air-sea interaction mechanism for the northward propagation of BSISO over the WNP (e.g., Roxy and Tanimoto 2012; Wang et al. 2018; Gao et al. 2019). Therefore, we calculated the bimodal ISO index (Kikuchi et al. 2012; Kikuchi 2020) and examined BSISO phases during the strong MT. Figure 16 is the diagram demonstrating the region where the BSISO convection is active and shows the region at lag = −5 – +5. In observational data (Fig. 16a), the BSISO convection is located over the Maritime continent (phase 5) at lag = −5 (red marker), and then moves to the WNP (phase 6) at lag = +5 (blue marker), demonstrating that convections associated with the active MT are related to the BSISO convection. Similarly, most members in each numerical experiment show that BSISO convections are located over the WNP when the MT becomes active (Figs. 16b–e). These results suggest that BSISO may affect the fluctuation of the MT, which is consistent with the findings of Zhou and Chan (2005).

BSISO phases (Kikuchi et al. 2012; Kikuchi 2020) during the strong MT for (a) NOAA, (b) NICOCO, (c) NICAM-A (ERA5), (d) NICAM-A (NOAA), and (e) NICAM-N. The texts and numbers in (a) represent the location and phase of major BSISO convections, respectively. The red (blue) marker in each figure indicates lag = −5 (+5).
There are many previous studies that have examined other mechanisms for the northward propagation of BSISO. Although this study does not focus on BSISO, similar mechanisms likely influence the northward propagation of convections following the MT, as BSISO is related to the MT. Jiang et al. (2004) proposed that mean easterly VWS and moisture advection in the PBL facilitate northward propagation of convections, with the latter being related to mean southerly wind and the meridional gradient of specific humidity in the PBL. In addition, Bellon and Sobel (2008) emphasized that mean northerly VWS also contributes to the northward propagation, while Dixit and Srinivasan (2011) indicated that mean zonal (meridional) VWS plays a role in the direction (speed) of the northward propagation. Given these previous studies, we simply compared the background field between NICOCO and NICAM-A (ERA5) (Fig. 17). Both NICOCO and NICAM-A (ERA5) realistically simulate mean VWS (Figs. 8b, c). Figure 17a shows easterly VWS around Philippines is larger in NICAM-A (ERA5) than in NICOCO. Figures 17b and 17c demonstrate that differences of vertical meridional wind shear and southerly in the PBL around Philippines are small. Figure 17d shows that the meridional gradient of specific humidity in the PBL is steeper in NICAM-A (ERA5) than in NICOCO while the gradients in NICAM-A (NOAA) and NICAM-N are nearly same as in NICOCO. Thus, the fact that NICAM-A (ERA5) simulates the meridionally steeper gradient of specific humidity may be attributed to the underestimation of SST (Fig. 6). Based on these results, the background is more favorable for the northward propagation of convections in NICAM-A (ERA5) than in NICOCO, although convections following the MT propagate faster in NICOCO, which further supports the contribution of air-sea interaction to the speed of the northward propagation. Comparing the result between NICOCO and NICAM-A (NOAA), NICOCO simulates slightly steeper specific humidity gradient than NICAM-A (NOAA), but the difference is smaller compared to the difference between NICOCO and NICAM-A (ERA5). Therefore, specific humidity in the PBL may be more contributable in NICOCO than in NICAM-A (NOAA), although the contribution is not expected to be essentially important to the difference of the MT fluctuation.

Comparisons of background fields. (a–c) show the difference between NICOCO and NICAM-A (ERA5) for (a) mean vertical zonal wind shear, (b) mean vertical meridional wind shear, and (c) mean southerly wind in the PBL. (d) shows mean specific humidity in the PBL averaged from 110°E to 130°E for ERA5 and simulations results. Hatched regions in (a–c) indicate that the difference is statistically significant at the 95 % confidence level based on Welch’s t-test. White regions in (a, b) and (c) indicate terrain at 850 hPa and 900 hPa, respectively.
In the NICAM-N experiments, while the northward propagation speed of the MT is similar to results in the NICOCO experiments (Fig. 13), surface wind and latent heat flux are slightly larger, and net surface heat flux is more negative than in NICOCO (Fig. 15, Table 4). This may result from the difference between air-sea coupled and uncoupled processes and the non-linearity of the related convection events. Assuming that SST conditions between NICOCO and NICAM-N are nearly same, heat flux continues to be provided to the atmosphere in NICAM-N since the prescribed SST does not change simultaneously even if the atmosphere attempts to influence the ocean, which suggests that air-sea coupling is important not to overestimate the intensity of convections. In addition to the coupling effect, the atmospheric internal variability would be closely related to the difference of heat flux between NICOCO and NICAM-N. For instance, because the BSISO phases are not always in the same phases as seen in Fig. 16, SST anomalies and their effects during the active MT may differ even if the same SST is prescribed. In other words, if positive (negative) SST anomaly covers the active MT region, much (less) heat flux would be provided to the atmosphere with much (less) convection. Because the internal variabilities may be canceled by ensemble mean, the response is likely to exhibit a linear relationship with SST as the number of ensemble members increases.
We investigated the role of air-sea interaction in the intra-seasonal fluctuation of the MT by comparing a global high resolved numerical model “NICAM” and its ocean coupled version “NICOCO.” We conducted 10-member ensemble numerical experiments with NICAM and NICOCO from 10 July to 1 September 2020. In the NICOCO simulations, the ocean was initialized with ocean reanalysis data, and a flux adjustment method was applied to achieve a realistic ocean state. The NICAM experiments were forced by two types of SST: one fluctuates only along the climatological seasonal cycle (SST-A), and the other is obtained from the results simulated in NICOCO (SST-N). NICAM experiment forced by SST-A derived from ERA5 was named as NICAM-A (ERA5), SST-A derived from NOAA as NICAM-A (NOAA), and SST-N as NICAM-N.
NICOCO successfully simulated a realistic ocean field; all the experiments represent small eastward extension of the MT, characteristic under a La Niña condition. However, the westerly and convections associated with the MT were more active in the NICAM-N experiments. Based on composite analysis during the strong MT, it is indicated that NICAM-A significantly underestimates the northward propagation speed of the westerly region and convections following the MT compared to NICOCO. This discrepancy can be attributed to differences of SST fluctuations. Over the WNP, when convective activity becomes active and many clouds are generated, the downward shortwave radiation is reflected by the clouds, reducing ocean surface heating. In addition, cyclonic circulations following convections accelerate (decelerate) the background westerly/southwesterly on the southern (northern) portion of the convections, which induces greater latent heat flux release on the southern side. Consequently, the net downward flux at the surface is more negative on the southern portion than on the northern portion of convections, leading to a dominant negative SST anomaly in the south and creating a meridional asymmetric SST anomaly structure. This SST fluctuation is not represented in NICAM-A and is inconsistent with surface flux fluctuations since SST is prescribed (i.e., such an SST fluctuation is not considered). As discussed in Section 4, the background field tends to be more favorable for the northward propagation of convections in NICAM-A (ERA5) than in NICOCO, although the propagation speed is faster in NICOCO than in NICAM-A (ERA5). Our results suggest that the difference in SST anomaly structure induced by the MT may result in the difference in the northward propagation speed. Therefore, an atmosphere-only forecast simulation would not represent the realistic northward propagation of the MT because SST fluctuations due to the atmosphere are not considered. The northward propagation can be successfully simulated if SST data influenced by the MT are used, although the intensity of convections is overestimated.
There are some issues to be addressed for conducting numerical experiments with NICAM and NICOCO. First, a flux adjustment (Masunaga et al. 2023) is necessary for simulating a realistic ocean field because there is SST drift influenced by the underestimation of cloud cover in NICAM (Kodama et al. 2021). The bias in NICAM needs to be reduced to achieve realistic simulations without flux adjustment. Next, we examined the fluctuation of the MT with high-resolution NICAM and NICOCO experiments in this study. However, the horizontal resolution of the atmosphere is 14 km, and this resolution may still be too coarse to investigate convections associated with the MT. With recent advancements in supercomputers, it has become feasible to conduct global numerical experiments at km-scale resolution (Stevens et al. 2019). Conducting such simulations will expand our understanding of the MT. Additionally, further analysis is required to fully comprehend the role of airsea interaction in the MT fluctuations. For instance, analyzing the moist static energy budget (e.g., Gao et al. 2019), the effect of the ocean mixed layer (e.g., Wolding et al. 2023), and the impact of ocean upwelling induced by tropical cyclones on SST fluctuations will provide deeper insights.
ERA5 data were downloaded from the Climate Data Store (https://doi.org/10.24381/cds.adbb2d47 and https://doi.org/10.24381/cds.bd0915c6). GLORYS12v1 data were downloaded from the Copernicus Marine Environment Monitoring Service (https://doi.org/10.48670/moi-00021). OLR and SST data were also downloaded from the NOAA (https://doi.org/10.7289/V5SJ1HH2 and https://doi.org/10.25921/RE9P-PT57).
Fig. S1. Time series of horizontal wind at 850 hPa averaged over 70–100°E for (a) ERA5, (b) NICOCO, (c) NICAM-A (ERA5), (d) NICAM-A (NOAA), and (e) NICAMN.
Fig. S2. Same as Fig. 13 except that the maximum of the westerly at 850 hPa over a narrower region [115–125°E, 10–15°N] is used as the composite metric.
Fig. S3. Same as Fig. 13 except that the minimum of OLR around the Philippines [110–130°E, 5–20°N] is used as the composite metric.
Fig. S4. The 95 % confidence range based on t-distribution for each variable in Fig. 14.
This research was supported by the JSPS KAKENHI Grants 22H01297/23K22568, 24K00707, and 24H02228, and the Ministry of Education, Culture, Sports, Science and Technology as “Program for promoting researches on the supercomputer Fugaku” JPMXP 1020200305, Large Ensemble Atmospheric and Environmental Prediction for Disaster Prevention and Mitigation. The ensemble simulations were performed on the supercomputer Fugaku (proposal numbers hp220167, hp230108, hp240179). Ocean spin-up simulations were performed on the supercomputer Wisteria/BDEC-01, provided by the Information Technology Center, The University of Tokyo (group ID gi55). ALERA3 data were provided by Dr. Akira Yamazaki. Constructive comments and suggestions from the editor and two anonymous reviewers are sincerely appreciated.