The role of air–sea coupling on November–April intraseasonal rainfall variability over the South Pacific

We investigate the impact of resolving air-sea interaction on the simulation of the intraseasonal rainfall variability over the South Pacific using the ECHAM5 atmospheric general circulation model coupled with the Snow-Ice-Thermocline (SIT) ocean model. We compare the fully coupled simulation with two uncoupled ECHAM5 simulations, one forced with sea surface temperature (SST) climatology and one forced with daily SST from the coupled model. The intraseasonal rainfall variability over the South Pacific is reduced by 17% in the uncoupled model forced with SST climatology and increased by 8% in the uncoupled simulation forced with daily SST, suggesting the role of air–sea coupling and SST variability. The coupled model best simulates the key characteristics of the two dominant patterns (modes) of intraseasonal rainfall variability over the South Pacific with reasonable propagation and correct periodicity. The spatial structure of the two rainfall modes in all three simulations is very similar, suggesting the dynamics of the atmosphere primarily generate these modes. The southeastward propagation of rainfall anomalies associated with two leading rainfall modes in the South Pacific depends upon the eastward propagating Madden–Julian Oscillation (MJO) signals from the Indian Ocean and western Pacific. Air-sea interaction improves such propagation as both eastward and southeastward propagations are substantially reduced in the uncoupled model forced with SST climatology. The simulation of both eastward and southeastward propagations considerably improved in the uncoupled model forced with daily SST; however, the periodicity differs from the coupled model. Such discrepancy in the periodicity is attributed to the changes in the SST-rainfall relationship with weaker correlations and the nearly in-phase relationship, attributed to enhanced positive latent heat flux feedbacks.


Introduction
The rainfall variability in the South Pacific shows a wide range of timescales: diurnal variability linked to radiative heating, local forcing, land-sea breeze, and large-scale circulations (Yang and Smith 2006;Kikuchi and Wang 2008;Cronin et al. 2015;Vincent and Lane 2016); synoptic variability associated with extratropical wave trains (Matthews 2012); Madden-Julian Oscillation (MJO) related intraseasonal variability (Matthews and Li 2005;Matthews 2012;Pariyar et al. 2019), and interannual and interdecadal variability influenced by El Niño-Southern Oscillation (ENSO) and Interdecadal Pacific Oscillation (IPO) (Folland et al. 2002;Vincent et al. 2011;Lintner et al. 2019). Among these, the intraseasonal timescale dominates the rainfall variability, explaining 50-60% (40-50%) of total daily rainfall variance along the South Pacific Convergence Zone (SPCZ) in November-April (May-October) (Pariyar et al. 2019). Two large-scale rainfall modes primarily dominate the intraseasonal rainfall variability in the South Pacific in both seasons (Matthews 2012;Pariyar et al. 2019). However, these leading modes are more pronounced in November-April, contributing up to 60% of daily intraseasonal rainfall variability over the large areas of SPCZ (Pariyar et al. 2019).
Our previous observational study identified two leading modes of intraseasonal rainfall variability in November-April over the South Pacific: the enhanced/suppressed and shifted SPCZ modes (Pariyar et al. 2019). These two modes have a typical timescale of 30-80-day, with dry and wet phases linked to the enhanced and suppressed phases of the Madden Julian Oscillation (MJO). They are connected to the MJO through two different mechanisms for southeastward propagation: wind-evaporation-SST and low-level frictional moisture convergence (see Supplementary Fig. 1). Interestingly, SSTs seem relevant in both cases, with positive SST anomalies typically leading the convection patterns. However, to what extent the intraseasonal SST anomalies interact with the atmosphere and contribute to the propagation of rainfall anomalies in the South Pacific is unclear.
One way to address this issue is to use general circulation models (GCMs). Many previous studies have discussed the role of air-sea interaction for the eastward propagation of intraseasonally varying rainfall anomalies, with a particular focus on the MJO in the Indian Ocean and western Pacific (Flatau et al. 1997;Hendon 2000;Woolnough et al. 2000;Kemball-Cook et al. 2002;Fu et al. 2003;Fu and Wang 2004;Sperber et al. 2005;Lin et al. 2006;Zhang et al. 2006;Pegion and Kirtman 2008;Wu et al. 2008;Bollasina and Nigam 2009;Wang and Seo 2009;DeMott et al. 2011;Roxy et al. 2013). Many modelling studies have shown that the MJO-induced SST variability over the tropical Pacific enhances low-level convergence, which increases the development of shallow convection and preconditions eastward propagation (Fu et al. 2003;Roxy et al. 2013;Tseng et al. 2015). Other similar studies argued that the local windevaporation-SST feedback increases surface evaporation, eventually increasing low-level moisture convergence (e.g., Marshall et al. 2008). In contrast, more recent studies suggest that positive SST feedback enhances the propagation of rainfall by amplifying horizontal moisture advection (DeMott et al. 2019;Zhou and Murtugudde 2020). However, studies on air-sea interaction in intraseasonal rainfall variability over the South Pacific are scarce. Further, we are not aware of any modelling studies on the southeastward propagation of intraseasonal rainfall variability over this region and their connection to the MJO. Therefore, we investigate these issues in detail using a state-of-art model.
Here we use the ECHAM5-SIT model, as it can realistically simulate the key characteristics of the MJO , which is strongly linked to the modes of intraseasonal rainfall variability in the South Pacific (Pariyar et al. 2019). Despite substantial improvements in GCMs, most are considered poor in their MJO simulation skill, and very few can realistically simulate the MJO (Kim et al. 2008;Hung et al. 2013;Jiang et al. 2015;Ahn et al. 2017Ahn et al. , 2020. ECHAM5 coupled with the snow-ice thermocline (SIT) onecolumn ocean model (ECHAM5-SIT) is one of the few models that reasonably simulates the MJO in terms of amplitude, periodicity, propagation speed, and vertical structure Jiang et al. 2015). We primarily use outputs from three model experiments performed by Tseng et al. (2015). These coupled and uncoupled experiments were used to demonstrate the key role of ocean-atmosphere interaction involving the upper few meters of the ocean in simulating the MJO. Here we use them to investigate how the representation of air-sea interaction and the simulation of the MJO impacts the intraseasonal variability over the South Pacific. We focus on November-April as the intraseasonal variability in the South Pacific is more pronounced, and the influence of the MJO is greatest (Pariyar et al. 2019).

Model, data, and methods
We use the simulations from ECHAM5-SIT as described in Tseng et al. (2015). The atmospheric component, ECHAM5, is an atmospheric general circulation model developed at the Max Planck Institute for Meteorology (MPI) (Roeckner et al. 2003). The model has a horizontal resolution of T63 (~ 1.8°) and 31 vertical layers from surface to 10 hPa. The Nordeng convection scheme was used in the model (Nordeng 1994). The ocean component is the one column snow-ice-thermocline (SIT) ocean model (Tu and Tsuang 2005;Tsuang et al. 2009). SIT is based on the turbulent kinetic energy approach suggested by Gaspar et al. (1990) and determines the ocean temperature, salinity and zonal and meridional currents in each depth according to the salinity, momentum, and energy budget in a one-column model. The ocean model used in this study has already been validated over the tropical ocean (Tu and Tsuang 2005), the South China Sea (Lan et al. 2010), and the Caspian Sea (Tsuang et al. 2001;Arpe et al. 2018). SIT has 42 vertical layers, with 12 in the upper 10 m of the ocean in its standard configuration. The fine upper ocean resolution enables air-sea interaction on diurnal to intraseasonal time scales to be well represented. To account for neglected horizontal processes, the ocean was weakly nudged (with a 30-day time scale) to the observed climatological ocean temperature below the 10-m depth.
We mainly use three experiments from Tseng et al. (2015) to investigate the role of SST and the impact of air-sea coupling on the simulation of intraseasonal variability in the South Pacific (Table 1). The first experiment (C-CTL) includes a fully coupled control simulation for 25 years. Note that the coupling was done only over the tropical domain (30° S-30° N), and outside of this region, observed SST climatology is prescribed. The second experiment (A-CLIM) comprises an atmosphere-only simulation forced with monthly SST climatology from the coupled simulation. The comparison of A-CLIM with C-CTL assesses the impact of tropical air-sea coupling on intraseasonal variability. The third experiment (A-DAY) is an uncoupled simulation forced with daily SST from the coupled simulation. Here, the ocean is forcing the atmosphere, but the 1 3 atmospheric influence on the ocean is absent. Even though the atmospheric feedback to the ocean is suppressed, the SST still varies on intraseasonal time scales and allows us to understand the impact of SST variability on the simulation of the intraseasonal variability. For more details on different experiments, readers are referred to Tseng et al. (2015).
Apart from the three main experiments, we use five additional simulations. To assess the role of ocean vertical resolution, we use two experiments in which the fine resolution of the upper ocean was replaced with the first layer of either 16.8 m (C-17 m) or 59.3 m (C-59 m) thickness. Precisely, the ocean layers between 0.05 m and 10 m (43.6 m) were removed and replaced by a single layer with a thickness of 16.8 m (59.3 m) for C-17 m (C-59 m). To assess the role of air-sea coupling in different regions, we use three experiments with regional coupling over the Indian Ocean (C-IO; 30° N-30° S, 50° E-100° E), the western Pacific (C-PO; 30° N-30° S, 110° E-180° E), and the Indian-western Pacific region (C-IPO; 30° N-30° S, 40° E-180° E).
We use the daily TRMM 3B42 version 7 rainfall estimates (Huffman et al. 2007) to compare the mean state and the Spatio-temporal variability of dominant large-scale rainfall patterns with the control simulation. We obtain the TRMM data with a horizontal resolution of 0.25° from 1998 to 2018. To compare the propagation mechanism, we use daily ERA-Interim data sets from 1998 to 2018 (Dee et al. 2011) with a horizontal resolution of 0.75 degrees. We interpolate both TRMM and ERA-Interim data sets to a regular 1.8 degrees grid to be consistent with the model.
The dominant intraseasonal rainfall modes are characterized by an Empirical Orthogonal Function (EOF) analysis based on 10-90-day bandpass filtered data for November to April similar to our previous study (Pariyar et al. 2019). Our entire analysis uses November to April as this is the rainy season in the South Pacific. Before EOF analysis, we compute the daily rainfall anomalies by subtracting the annual rainfall climatology. We detrend the anomalies and apply a 10-90-day Butterworth bandpass filter to remove the longterm trend and retain intraseasonal variability. We apply a 30-80-day filter to the principal components (PCs) as it is a typical MJO time scale (Wheeler and Hendon 2004) as well as the dominant time scale for the two EOF modes in observations (Pariyar et al. 2019). We perform the EOF analysis over the western Pacific (120° E-150° W; 20° N-25° S) to be consistent with our previous study (Pariyar et al. 2019). Note that all the analysis presented for observations is similar to our previous study unless stated otherwise. Readers are referred to Pariyar et al. (2019) for more details. We repeat a similar analysis for the model data. To compare the observed EOF modes with the model simulations, we project the model simulations onto the observed EOF modes. The main reason for doing this is to ensure that the EOF modes in both model and observations represent the same mode of variability.
We test the significance of the leading EOF modes against a stochastic null hypothesis based on first-order autoregressive (AR1) and spatial isotropic diffusive processes (Dommenget 2007). This test compares the observed EOF modes with the EOF modes of the null hypothesis fitted to the observations. The observed EOF modes that are statistically different from the fitted EOF modes are considered distinct EOF modes. We also apply the North's rule to distinguish the leading modes from the higher modes (North et al. 1982).
We compute the power spectra for two principal components (PCs) for November-April and then average over 25 years. We test for statistical significance of the power spectra from a red-noise spectrum at the 95% confidence level (Torrence and Compo 1998). We compute the 95% confidence level by multiplying the red noise spectrum by the 95 th percentile value of the chi-squared distribution (Gilman et al. 1963). We compute the lag-correlation between two PCs to understand the phase relationship between leading EOFs. We apply the two-tailed Student's t-statistic at the 95% confidence level for the statistical significance of the correlation.
We compute lagged composites to represent the lifecycle and propagation characteristics of the leading intraseasonal rainfall modes. The composites are based on intraseasonal rainfall events identified by standardized PCs with values higher than 0.5 standard deviations. We average the data for days before and after the event to illustrate the evolution of rainfall events. We compute the lag-composites for other variables similarly. For the statistical significance of the lag-composites, we implement the two-tailed Student's t-statistic at the 95% confidence level (Harrison and Larkin 1998). To represent the eastward propagation over the Indian Ocean and the western Pacific, we average the rainfall composites from 5° N to 10° S; these composites represent the eastward propagation from the Indian Ocean to the western Pacific. We compute the southeastward propagation for the South Pacific by averaging rainfall anomalies along the diagonal axis corresponding to the SPCZ in each EOF mode. We construct the diagonal axis by identifying the maximum rainfall anomaly for each latitude along the SPCZ in each EOF mode. Later, we fit a linear diagonal line based on these maximum rainfall values. We average the data from 10° to the east to 10° to the west of the diagonal axis for each latitude.
To understand the phase relationship between leading EOF modes, we compute the lag-correlation coefficients between two PC time series. Likewise, we perform a lagcorrelation analysis to establish the phase relationship between SST and other atmospheric variables, namely latent heat flux (LHF, positive upwards), sensible heat flux (SHF, positive upwards), rainfall (PCP), 10-m wind speed (WS). We perform this analysis for both coupled and uncoupled simulations. To establish a relationship between MJO and the leading EOF modes, we compute the lag-correlation coefficients between the real-time multivariate MJO (RMM) index (Wheeler and Hendon 2004) and the PC time series. We obtain the daily RMM indices (RMM1 and RMM2) from the Australian Bureau of Meteorology website (http:// www. bom. gov. au/ clima te/ mjo/). We compute the RMM indices as described in Wheeler and Hendon (2004) for the model. We apply the two-tailed Student's t-statistic at the 95% confidence level for the statistical significance of lag correlation coefficients.
To quantify the LHF-rainfall relationship, we compute a feedback parameter for latent heat flux based on the linear regression between daily anomalies of latent heat flux and rainfall (Dellaripa and Maloney 2015;Bui et al. 2020). Note that both LHF and rainfall anomalies are expressed in W m −2 units. The feedback parameter quantifies the positive feedback between surface wind, LHF, and rainfall and is thus useful to quantify the role of LHF feedback in coupled and uncoupled models.

Mean and intraseasonal variability of the rainfall
The mean climate from coupled and uncoupled simulations with an identical model experiment has been evaluated previously by Tseng et al. (2015) over the Indo-Pacific region; therefore, we will only briefly summarize the comparison of the mean state and intraseasonal (10-90-day) variability of rainfall from the coupled simulation with observations over the South Pacific ( Fig. 1). ECHAM5-SIT simulates the similar spatial distribution of the mean rainfall over the South Pacific. The control simulation generally overestimates the rainfall climatology with a slightly more zonal orientation of the SPCZ (black diagonal dashed lines, Fig. 1b). Also, the component of the ITCZ to the north of the equator is largely absent. The more zonal orientation of the SPCZ is a common problem in most global climate  (Brown et al. 2011) and is linked to SST biases (Brown et al. 2011;Shen et al. 2016). The coupled model used in our study also has a warm SST bias over the South Pacific , and such a warm SST bias can explain these rainfall biases. Despite the mean state bias, the ratio of intraseasonal variance to the total variance along the SPCZ in the coupled model is similar to observations (Supplementary Fig. 2). Interestingly, the intraseasonal variance ratio in two uncoupled models is similar to that of the coupled model, suggesting that the model's mean state errors do not greatly affect its simulation of intraseasonal variability. Therefore, we believe that despite the mean state biases, ECHAM-SIT is a useful tool to study intraseasonal rainfall variability in the South Pacific, as it is for studying the MJO Chang et al. 2015;Jiang et al. 2015). C-CTL nicely captures the observed spatial pattern of intraseasonal variability (Fig. 1f); however, the overall intraseasonal variability is slightly weaker in the model. Past studies have shown a poor simulation of the tropical intraseasonal variability is linked to the simulation of the basicstate surface westerlies in the Indian Ocean and western Pacific (Hendon 2000;. The climatology of the 10-m zonal winds over the western Pacific in the coupled and uncoupled model simulations is similar to observations ( Supplementary Fig. 3). In particular, the westerly winds over the western Pacific are reasonably well simulated by the model. However, a key difference is that the westerly flow in the coupled model is slightly weaker and extends farther eastward than in observations. Nevertheless, the overall agreement with observations in terms of these features further supports the usefulness of ECHAM-SIT for this kind of study.

Dominant intraseasonal (10-90-day) rainfall modes
We now show that ECHAM5-SIT can capture the observed Spatio-temporal structure of intraseasonal rainfall variability in the South Pacific. We perform an EOF analysis of the daily 10-90-day filtered rainfall anomalies over the western Pacific (25° S-20° N to 120° E-150° W) for observations and coupled model. The first two EOF modes explain about 11% of the 10-90-day rainfall variability in observations, whereas, in the coupled model, this value is about 15%. The coupled model captures the large-scale spatial structures of the leading EOF modes reasonably well (Fig. 2d, e). The EOF1 mode in observations has maximum rainfall variance along the mean SPCZ position (Fig. 2a), whereas the EOF2 mode has maximum southwest of the climatological SPCZ position (Fig. 2b). We term the EOF1 (EOF2) mode as enhanced (shifted) SPCZ mode (Pariyar et al. 2019), and these features are well represented in the model. Despite these similarities, there is one main difference. In C-CTL, the major axis of rainfall anomalies in the SPCZ in both EOF modes is more zonal than observations (diagonal black lines, Fig. 2d, e), consistent with the more zonal orientation of the mean SPCZ (Fig. 1b). The comparison of power spectra for two EOF modes between C-CTL and observations shows comparable results with maximum variance concentrated in the 30-80-day period. The power spectra of two EOF modes in CCTL peak around 40-50 days (Fig. 2f), similar to observations (Fig. 2c). However, the variance in the PCs is lower in the coupled model, consistent with the weaker intraseasonal variability in the model. The two PCs in C-CTL show good correspondence in the 30-80-day period (r ~ 0.42; lag ~ 10 days) with PC1 leading PC2, suggesting that the two EOF modes may partly represent the same mode of variability, unlike observations (r ~ 0.23; lag ~ 8 days).
The MJO-South Pacific teleconnection in C-CTL agrees well with observations (figure not shown). In observations, RMM1 leads PC1 by two days (r ~ 0.56) and negative of PC1 leads RMM2 by eight days (r ~ − 0.58); and RMM1 leads PC2 by 12 days (r ~ 0.41), and PC2 leads RMM2 by one day (r ~ 0.46). These values are quite comparable to C-CTL, where RMM1 leads PC1 by four days (r ~ 0.66) and negative of PC1 leads RMM2 by six days (r ~ 0.49); and for EOF2, PC2 leads RMM1 by five days (r ~ 0.54), and RMM2 leads PC2 by three days (r ~ 0.39).

Rainfall propagation characteristics
We summarize the propagation features in time-longitude and time-latitude plots based on the lag-composites (Fig. 3). The time-longitude plots represent the eastward propagation from the Indian Ocean to the western Pacific. Similarly, the time-latitude plots represent the southeastward propagation in the South Pacific. We consider the eastward propagation over the Indo-Pacific region because the southeastward propagation of the rainfall anomalies on the intraseasonal time scale over the South Pacific is associated with the eastward propagating MJO signal from the Indian Ocean and the western Pacific (Pariyar et al. 2019).
The eastward propagation of the rainfall anomalies from the Indian Ocean to the western Pacific in both EOF modes simulated in C-CTL is comparable to observations (Fig. 3e,  f); however, the amplitudes of the rainfall anomalies are slightly weaker in the model, particularly over the Indian Ocean and western Pacific. The southeastward propagation over the South Pacific is also well simulated in C-CTL. In Fig. 3 Time-longitude and time-latitude composites of 30-80-day rainfall anomalies (mm day −1 ) for a-d observations, e-h coupled model, i-l uncoupled model forced with SST climatology, and m-p uncoupled model forced with daily SST. The first two columns represent the eastward propagation from the Indian Ocean to the western Pacific, and the last two columns represent the southeastward propagation in the South Pacific. For the eastward propagation, we compute the average between 5° N and 10° S. We average along the diagonal axis over the South Pacific (black diagonal dashed lines, Fig. 2) for the southeastward propagation. The diagonal axis is constructed by following the maximum rainfall anomalies along the rainfall bands over the South Pacific associated with the EOF modes in observations. For each latitude, the region from 10 degrees to the east to 10 degrees to the west of the diagonal axis is chosen for averaging, and the respective longitudes are also provided. The rainfall lag composites are computed from standardized PCs for selected strong events with values greater than 0.5 standard deviations. Only significant values at the 95% confidence level are shown ◂ 1 3 observations, the southeastward propagation between 10° N and 10° S is quite pronounced with relatively stationary signals between 10 and 20° S (Fig. 3c, d). This feature is nicely simulated in C-CTL with comparable amplitudes of rainfall anomalies (Fig. 3g, h).

Mean state and intraseasonal variability of the rainfall
The spatial structures of the daily rainfall climatology in both uncoupled simulations closely resemble the coupled simulation with comparable magnitudes (Fig. 1c, d). However, there is one noticeable difference. The orientation of the SPCZ in both uncoupled models is more zonal (diagonal dashed lines, Fig. 1). Consequently, the location of the maximum rainfall is shifted more to the equator. These differences in the mean SPCZ positions in both uncoupled models can be attributed to air-sea coupling, as they have identical SST climatology to C-CTL. The overall intraseasonal rainfall variability in the South Pacific in A-CLIM is reduced compared to C-CTL (Fig. 1g). In contrast, the intraseasonal variability is increased in A-DAY (Fig. 1h). In particular, the changes in intraseasonal variability are more evident over the SPCZ. For quantitative assessment, we compute weighted area average standard deviations over the SPCZ (160° E-170° W, 5° S-15° S) for all three simulations (black rectangular box, Fig. 1h). We chose this region because the maximum standard deviations are mainly concentrated along the SPCZ in all three simulations. There is a 17% reduction of the intraseasonal variability in A-CLIM compared to C-CTL. Such a reduction in the intraseasonal variability can be attributed to suppressed SST variability and the absence of air-sea coupling .
We observe an 8% increase in the intraseasonal rainfall variability in A-DAY compared to C-CTL, suggesting that the prescribed daily SSTs inflate rainfall variability compared to coupled simulation. Pegion and Kirtman (2008) observed a similar increase in rainfall variability in the uncoupled model forced with daily SST from a coupled model. They further argue that, in the uncoupled model, the prescribed SST can cause local rainfall anomalies through positive latent heat feedback if the wind circulation is favourable for convection. In contrast, in reality, and in the coupled model, the SST anomalies are damped by the atmospheric response. Thus, relatively larger rainfall anomalies are expected in the uncoupled simulation, so the rainfall variability increases.

Simulation of leading EOF modes and their relation to the MJO
The first two EOF modes in both uncoupled models explain roughly 14 and 15% of the total 10-90-day rainfall variability, comparable to the coupled model (15%). The spatial structures of the two leading EOF modes in both uncoupled simulations are similar to the coupled simulation ( Fig. 2g, h, j, k). However, the maximum rainfall Fig. 4 Time-longitude and time-latitude composites of 30-80-day rainfall anomalies (shaded, mm day −1 ) and other atmospheric fields [contours, sea surface temperature (°C, contour interval 0.05, beginning at ± 0.05); sea level pressure (hPa, contour interval 0.2, beginning at ± 0.2); 1000-hPa moisture flux divergence (10 -6 g kg −1 s −1 , contour interval 6, beginning at ± 5); latent heat flux (w m −2 , contour interval 4, beginning at ± 4, positive upwards); sensible heat flux (w m −2 , contour interval 0.5, beginning at ± 0.5, positive upwards); 10-m wind speed (m s −1 , contour interval 0.3, beginning at ± 0.3] for a-f coupled model, and g-l uncoupled model forced with daily SST. We average along the diagonal axis over the South Pacific (black diagonal dashed lines, Fig. 2). The diagonal axis is constructed by following the maximum rainfall variance along the rainfall bands over the South Pacific associated with each EOF mode in observation. For each latitude, the region from 10° E and 10° W of the diagonal axis is chosen for averaging, and the respective longitudes are also provided. The lag composites are computed from standardized PCs by selecting strong events with values greater than 0.5 standard deviations. Only significant values at the 95% confidence level are shown. Solid contours represent positive anomalies, and dashed contours are for negative anomalies Although the overall large-scale rainfall structures for both EOF modes in both uncoupled models are comparable to the coupled model, the power spectra differ. The power spectra of two PCs for A-CLIM show no evidence of pronounced peaks in the 30-80-day period (Fig. 2i). The overall spectra for both PCs are relatively flat without any preferred oscillatory period. In contrast, the power spectra of two PCs in A-DAY show pronounced peaks in the 30-80-day period with comparable amplitudes to C-CTL (Fig. 2l); however, the peaks with maximum power are shifted to the left from about 40 days in C-CTL to 60 days in A-DAY.
To further quantify the impact of air-sea coupling on the variance associated with two PCs, we compute the mean power associated with two PCs by averaging the power spectrum in the 30-80-day band (Pow SP ). Note that each PC's mean power is computed separately and averaged to obtain a single value. In the absence of air-sea coupling and SST variability, the variance associated with two PCs in the 30-80-day band is significantly reduced in A-CLIM (Pow SP ~ 0.19) than C-CTL (Pow SP ~ 0.34). Interestingly, the variance in A-DAY (Pow SP ~ 0.33) is quite comparable to C-CTL.
The MJO-South Pacific teleconnection also differs markedly in A-CLIM compared to C-CTL, and A-DAY, based on lag-correlation between PCs and RMM indices (figure not shown). The maximum MJO-PCs lag-correlation coefficients in the 30-80-day period in C-CTL range from -0.39 to + 0.66. These values are reduced to − 0.26 to + 0.3 in A-CLIM and somewhat comparable in A-DAY (r ~ − 0.46 to + 0.7), consistent with the better simulation of the MJO in A-DAY compared to A-CLIM .

Propagation characteristics and mechanism
The 30-80-day rainfall propagation characteristics in the uncoupled simulations are summarized in time-longitude and time-latitude plots (Fig. 3i-p). The eastward propagation from the Indian Ocean to the western Pacific is significantly reduced in A-CLIM (Fig. 3i, j). Consistent with this, the southeastward propagation over the South Pacific is also considerably weakened (Fig. 3k, l), attributed to the degraded simulation of the MJO in the absence of air-sea coupling and SST variability  and its teleconnection to the South Pacific, as discussed in Sect. 4.2.
In A-DAY, both eastward and southeastward propagations are improved (Fig. 3m-p). In particular, the eastward propagations from 110° E to 150° W and southeastward propagations from 10° N to 20° S are comparable to C-CTL for both EOF modes. Such improvements in both eastward and southeastward propagations suggest the role of intraseasonal SST anomalies for rainfall propagation over the Maritime Continent and the South Pacific, partly attributed to relatively better representation of the MJO in A-DAY compared to A-CLIM .
The observed southeastward propagation mechanisms associated with the two EOF modes are wind-evaporation-SST feedback and low-level frictional moisture convergence for EOF1 and EOF2 modes (see Supplementary  Fig. 1 and Pariyar et al. 2019 for more details). However, the coupled model simulates a low-level frictional moisture convergence propagation mechanism for both EOF modes that resemble EOF2 in observations (Fig. 4). In particular, warm SST anomalies lead maximum rainfall by 12-14 days (Fig. 4a), and low pressure and low-level moisture flux convergence anomalies slightly lead the rainfall (Fig. 4b, c). The evaporation is enhanced as represented by positive latent heat flux anomalies (Fig. 4d), and surface wind speeds are increased 2-3 days after the maximum rainfall (Fig. 4f). Unlike in observations, the wind-evaporation-feedback seems less relevant for the southeastward propagation of rainfall anomalies associated with the EOF1 mode as maximum wind speed and positive anomalies of surface fluxes are either in phase or follow rainfall for both EOF modes in the model.
The relationships between rainfall and other atmospheric variables in A-DAY are similar to C-CTL ( Fig. 4g-l); however, the SST-rainfall relationship differs. In particular, the warm SST anomalies are more or less in phase with rainfall anomalies (Fig. 4g). The shorter lag between SST and rainfall on the intraseasonal timescale is a key characteristic of uncoupled models (Fu et al. 2003;Fu and Wang 2004;Zheng et al. 2004;Pegion and Kirtman, 2008;DeMott et al. 2014). The surface fluxes modulated by the atmosphere have no impact on the SST in the uncoupled model. Therefore the atmospheric adjustment is almost in phase with the prescribed SST anomalies. The correct SST-rainfall phase relationship is essential because the warm SST ahead of the convection acts to destabilize the atmospheric boundary layer and promote frictional convergence (Flatau et al. 1997;Waliser et al. 1999;Kemball-Cook et al. 2002;Maloney and Kiehl 2002;Benedict and Randall 2011;Tseng Fig. 6 Time-longitude and time-latitude composites of 30-80-day rainfall anomalies (mm day −1 ) for a-d C-17 m, e-h C-59 m, i-l C-IO, m-p C-IPO, and q-t C-PO. For the eastward propagation, we compute the average between 5° N-10° S. We average along the diagonal axis over the South Pacific (black diagonal dashed lines, Fig. 2) for the southeastward propagation. The diagonal axis is constructed by following the maximum rainfall anomalies along the rainfall bands over the South Pacific associated with each EOF mode in observations. For each latitude, the region from 10° E and 10° W of the diagonal axis is chosen for averaging, and the respective longitudes are also provided. The rainfall lag composites are computed from standardized PCs by selecting strong events with values greater than 0.5 standard deviations. Only significant values at the 95% confidence level are shown ◂ et al. 2015). Therefore, it is interesting to investigate how the amplitude of the SST-rainfall relationship changes between coupled and uncoupled simulations and how such changes are attributed to changes in the SST relationship with the surface fluxes and wind speed.
We perform a lag-correlation between SST and other atmospheric variables (PCP, LHF, SHF, WS) (Fig. 5). We estimate the lag correlation coefficients for C-CTL and A-DAY by taking a weighted area average over 160° E-170° W and 5° S-15° S (black rectangular box, Fig. 1h). We choose this region for two reasons. First, it corresponds to the region of maximum intraseasonal variance. Second, it covers the region of maximum variance corresponding to two EOF modes in both coupled and uncoupled simulations (Fig. 2d, j).
The phase relationship between SST and rainfall has changed substantially in A-DAY compared to C-CTL (Fig. 5a). The correlation coefficients between SST and rainfall are 0.66 and 0.44 for C-CTL and A-DAY, respectively, with positive lags of 13 and 3 days. The weaker SSTrainfall relationship in A-DAY is supported by a weaker association between SST and other variables (SHF, LHF, WS). We observe the highest reduction in the maximum positive correlation coefficients between SST and wind speed (Fig. 5d). The weaker and nearly in phase SST-wind speed and SST-LHF/SHF relationships (Fig. 5b, c) indicate a lack of atmospheric influence on the ocean as the perturbation of SST anomalies over this region of the tropical Pacific is largely forced by wind-driven surface latent heat flux anomalies (Wu and Chen 2015). In the absence of atmospheric feedback to local SST, the rainfall anomalies in A-DAY would grow as long as the large-scale wind anomalies and associated surface heat fluxes are favourable, as indicated by significantly higher LHF feedback parameter (Fig. 5e), based on linear regression of LHF and rainfall anomalies (Dellaripa and Maloney 2015;Bui et al. 2020). The LHF anomalies at maximum positive lag are roughly 6% of the rainfall anomalies for C-CTL, which is largely exaggerated in A-DAY with a value of 9%. This is why the intraseasonal variability over the SPCZ is higher in A-DAY than in C-CTL (Fig. 1f, h).

Role of ocean vertical resolution and regional coupling
The analysis above shows the importance of air-sea interaction for the propagation and timescale of the intraseasonal variability in the South Pacific. To further assess the role of local versus remote air-sea interaction and the role of the upper ocean vertical resolution, we use five additional experiments: two experiments with coarser ocean vertical resolution (C-17 m and C-59 m) and three regional coupling simulations (C-IO, C-IPO, and C-PO). Both experiments with coarser ocean vertical resolution, in general, capture both eastward and southeastward propagations of rainfall anomalies and are somewhat comparable to C-CTL ( Fig. 6a-d, e-h). Interestingly, the eastward propagation in C-59 m is slightly slower than in C-17 m, and the southeastward propagation over the South Pacific is relatively weak. For regional coupling experiments, both eastward and southeastward propagation of rainfall anomalies is degraded in all three experiments (Fig. 6i-t). Among three simulations, C-IPO best simulates both eastward and southeastward propagations (Fig. 6m-p). To compare the performance of eight simulations quantitatively, we use six metrics: three for the MJO, two for the intraseasonal variability over the South Pacific, and one for MJO-South Pacific teleconnection. The first metric is the eastward/westward power ratio (EWP MJO ), obtained by dividing the sum of spectral power over the 30-80-day band by its westward propagating counterpart. This metric estimates the eastward propagating feature of the MJO (Kim et al. 2009). The second metric, a measure of the MJO periodicity (P MJO ), is formulated by computing twice the time interval between maximum and minimum lag-correlation between 30 and 80-day filtered RMM indices. The third metric represents the mean power in the 30-80-day band as represented by the power spectra of two RMM indices (Pow MJO ). The power spectrum is averaged over the 30-80day band for each RMM index and later averaged for both RMM indices. The two metrics representing intraseasonal variability over the South Pacific are the periodicity (P SP ) and mean power in the 30-80-day band (Pow SP ), similar to MJO indices. To define a metric that measures the teleconnection between the MJO and the two EOF modes in the 30-80-day period (Rmax MJO-SP ), we average the absolute values of maximum and minimum lead-lag correlation coefficients between PCs and two RMM indices. The correlation coefficients between each PC and two RMM indices are averaged to produce a single metric. The summary of metrics for observations and model experiments is presented in Table 2.
Among eight simulations, C-CTL best simulates the eastward propagating characteristics of the MJO with a value of the EWP MJO that is comparable to observations. The MJO propagation is somewhat captured in C-17 m, C-59 m, and C-IPO but degraded in uncoupled and regionally coupled experiments. Interestingly, the observed periodicity of the MJO in the 30-80-day period is well reproduced in all simulations except A-DAY. In A-DAY, the MJO periodicity is a little longer (P MJO ~ 46 days). The mean power in the 30-80-day is consistently higher than observations in all simulations; however, the value differs among models. For example, the mean power is largely reduced in A-CLIM and C-IO compared to C-CTL. The mean power is quite comparable to C-CTL for the remaining seven experiments. Even though the MJO periodicity and the mean power in the 30-80-day band in the majority of the experiments are similar, the peaks with maximum power for both RMM indices are more pronounced in C-CTL, A-DAY, C-17 m, C-59 m, and C-IPO (figure not shown).
The periodicity of the EOF modes in the 30-80-day period (P SP ) is pretty similar among all simulations except A-DAY. In A-DAY, the periodicity is a little longer (P SP ~ 50 days), consistent with the longer periodicity for the MJO. The observed mean power in the 30-80-day band (Pow SP ~ 0.45) is best simulated in C-CTL (Pow SP ~ 0.34); however, the mean power is largely reduced in all other experiments except A-DAY. Although the mean power in the 30-80-day band is quite similar among simulations, the peaks associated with the power spectra of two PCs are more pronounced in C-CTL, A-DAY, C-17 m, C-59 m, and C-IPO (figure not shown). It is noteworthy that the MJO related variability in the 30-80-day period not necessarily dictates the corresponding variance associated with two EOF modes in the South Pacific but rather enhances the oscillation of these modes. It is evident from Table 2 that the mean power in the 30-80-day band for the MJO is quite similar in all simulations except A-CLIM and C-IO; however, in the case of two EOF modes in the South Pacific, the mean power is comparable among models except C-CTL and A-DAY. Interestingly, the power spectra for two PCs in C-CTL, A-DAY, C-17 m, C-59 m, and C-IPO show pronounced peaks, consistent with the MJO.
Even though the characteristics of EOF modes, the phase relationship between PCs, and periodicity are comparable among simulations, their relationship with the MJO differs. For example, C-CTL best simulates the MJO-South Pacific teleconnection as indicated by comparable values of Rmax MJO-SP to observations. The coupled models, in general, better represent the MJO-South Pacific teleconnection than the uncoupled model. In particular, A-CLIM and C-IO perform worst in representing MJO-South Pacific teleconnection. The considerably weaker MJO-South Pacific teleconnection in A-CLIM, C-IO and C-PO is consistent with the degraded simulation of the MJO in these models.

Summary and discussion
We investigate the role of air-sea coupling and SST variability in the simulation of dominant intraseasonal (10-90-day) rainfall modes over the South Pacific using the ECHAM5-SIT model. The fully coupled ECHAM5-SIT model is capable of capturing the key characteristics of the intraseasonal rainfall variability over the South Pacific and its relationship with the MJO. It is, therefore, a useful tool for studying the mechanisms for the intraseasonal rainfall variability over this region.
Air-sea interaction improves the intraseasonal variability over the South Pacific. The overall intraseasonal rainfall variability along the SPCZ is reduced with suppressed air-sea interaction. The intraseasonal variability in the 30-80-day period is also reduced in the leading EOF modes over the South Pacific with suppressed air-sea interaction. However, a considerable fraction of intraseasonal variability is observed without air-sea interaction, attributed to the internal dynamics of the atmosphere.
The spatial structure of the dominant patterns of intraseasonal rainfall variability appears independent of air-sea interaction. The rainfall patterns associated with the two EOF modes are present in both coupled and uncoupled simulations and are comparable to observations. The spatial structure of the dominant modes remains nearly identical in coupled and uncoupled models, while their amplitudes differ. This strongly suggests that the internal dynamics of the atmosphere primarily generate these intraseasonal modes.
Prescribing intraseasonal SST variability improves both eastward and southeastward propagations associated with two EOF modes. In particular, daily varying SSTs produce rainfall variability associated with two EOF modes similar to the coupled model, but the period is longer. Such a longer periodicity in the absence of air-sea interaction is consistent with a longer periodicity of the MJO . The SST-rainfall relationship dramatically changes when the atmosphere does not interact with the ocean, leading to asymmetry in the SST-rainfall lead-lag relationship and longer oscillation of rainfall events.
Air-sea interaction and intraseasonal SST variability lead to the differences in the simulation of the MJO related eastward propagating rainfall variability, the southeastward propagation of rainfall anomalies over the South Pacific, and their teleconnection. In the absence of air-sea interaction, ECHAM5 fails to simulate the southeastward propagation over the South Pacific, consistent with significantly weaker eastward propagation over the Indian Ocean and western Pacific and degraded simulation of the MJO. On the other hand, prescribing SST variability improves both eastward and southeastward propagations. Ocean vertical resolution has a considerable impact on the simulation of these propagating features in the model as both eastward and southeastward propagations are degraded in the model with coarser vertical ocean resolution. Likewise, coupling over the entire Indo-Pacific basin improves both eastward and southeast propagations compared to locally coupled models, suggesting the importance of air-sea coupling in the whole Indo-Pacific basin.
The differences in southeastward propagation of rainfall anomalies over the South Pacific between coupled and uncoupled simulations are partly attributed to the eastward propagating signals in the Indian Ocean and western Pacific. These eastward propagating rainfall signals are primarily associated with the MJO and play an important role in enhancing the southeastward propagations in the South Pacific. Models with better MJO and associated eastward propagation tend to have more realistic southeastward propagation in the South Pacific and a better MJO-South Pacific teleconnection than models with degraded MJO simulation.
The southeastward propagation for two EOF modes in ECHAM-SIT appears to be associated with low-level moisture convergence that leads convection by 2-3 days.
However, in observations, the two EOF modes show different mechanisms, i.e., wind-SST-evaporation for EOF1 mode and low-level frictional moisture convergence mechanism for EOF2 mode (see Pariyar et al. 2019 and Supplementary Fig. 1). Thus, the observed wind-SST-evaporation mechanism is less relevant in the model. The simulation of propagation mechanisms for the tropical intraseasonal variability, particularly the MJO, is challenging, and the propagation mechanism differs among models (DeMott et al. 2015). Tseng et al. (2015), in their analysis with identical model experiments, argued that the improved MJO propagation in the coupled models is mainly due to enhanced frictional convergence ahead of the convection. In contrast, DeMott et al. 2019, based on four general circulation models, suggested that the improved MJO propagation in coupled experiments is mostly contributed by enhanced horizontal moisture advection. Therefore, it is interesting to investigate the role of moisture convergence and horizontal advection in more detail through a moist static energy budget. Likewise, similar experiments with other models are necessary to understand whether the propagation mechanisms simulated in ECHAM-SIT are model-dependent and why the wind-evaporation-feedback mechanism is not apparent in this model.
Our previous observational study hypothesized that the connection between the South Pacific and South China Sea for the EOF2 mode might explain the differences in the mechanisms for the two EOF modes (Pariyar et al. 2019). We observe statistically significant low-level wind circulation anomalies along the South China Sea for EOF2 mode, which is absent in the EOF1 mode. Therefore, we proposed that the large-scale forcing from the South China Sea could contribute to the low-level frictional convergence over the South Pacific. However, we do not find this connection in the coupled model as the circulation anomalies over the South China Sea are missing in both EOF modes ( Supplementary Fig. 4). Therefore, the difference in observed propagation mechanisms is not necessarily linked to the large-scale forcing from the South China Sea, but rather the absence of such a feature in our model could partly explain why the two modes are better related to each other and the MJO.
In summary, all our model simulations produce intraseasonal variability, but not all of them can produce southeastward propagating intraseasonal variability in the South Pacific. These differences can be attributed to air-sea coupling, intraseasonal SST variability, and representation of eastward propagating MJO variability in ECHAM-SIT. These findings motivate future work on intraseasonal variability in the South Pacific, particularly a detailed modelling study of propagation mechanisms based on moisture budget analysis.