The effect of diurnal warming of sea‐surface temperatures on the propagation speed of the Madden–Julian oscillation

The diurnal warm layer in the upper ocean develops during low surface winds and high incoming solar radiation conditions, often increasing sea‐surface temperatures (SSTs) by up to 1°C. The suppressed phase of the Madden–Julian Oscillation (MJO) favours the formation of such a layer. Here, we analyse the coupled ocean–atmosphere and atmosphere‐only numerical weather prediction systems of the UK Met Office to reveal that important differences arise from the representation of the diurnal warm layer in the coupled model. Though both models are skilful in predicting the MJO to at least a 7‐day lead time, the coupled model predicts approximately 12% faster MJO RMM phase speed propagation than the atmosphere‐only model due to the ability to resolve diurnal warming in the upper ocean that rectifies onto MJO‐associated SST anomalies. The diurnal warming of SST (dSST) in the coupled model leads to an increase in daily mean SST compared with the atmosphere‐only model persisted foundation SST. The strength of the dSST in the coupled model is modulated by MJO conditions. During suppressed MJO conditions on lead day 1, the dSST is enhanced, leading to 0.2°C warmer daily mean MJO‐associated SST anomalies and increased convection in the coupled model by lead day 7. During active MJO convection, the dSST is suppressed, leading to 0.1°C colder MJO‐associated SST anomalies in the coupled model and reduced convection by lead day 7. This variability in dSST further amplifies the MJO propagation speed, underlining the importance of the two‐way feedback between the MJO and the diurnal cycle of SST and the need to accurately represent this process in coupled models.


INTRODUCTION
The Madden-Julian Oscillation (MJO) is the main mode of intraseasonal (30-90 days) variability in the tropics (Madden & Julian, 1971).It is comprised of regions of enhanced and suppressed convection, O(10,000 km), propagating eastward with a phase speed of 5 m⋅s −1 across the tropics.The MJO convective anomalies typically originate in the west equatorial Indian Ocean (EIO), cross the Maritime Continent (MC; Indonesia, Philippines, and Papua New Guinea), and dissipate over the Pacific Ocean.
The MJO interacts with many global weather and climate patterns-for a review, see Zhang (2013)-and remains a challenge in subseasonal model predictability (e.g., Kim et al., 2019;Vitart, 2017).In the past decade, the rise of coupled ocean-atmosphere numerical weather prediction (NWP) and climate models has led to improvements in MJO predictions and simulations compared with atmosphere-only models (e.g., Ahn et al., 2017;Vitart, 2017).Current coupled ocean-atmosphere seasonal to subseasonal models predict the MJO out to 3-4.5 weeks (Kim et al., 2019); however, many underestimate its propagation speed and amplitude, especially over the MC (Kim et al., 2014;2019;Vitart, 2017;Xiang et al., 2015).Ocean-atmosphere feedbacks play an important role in MJO propagation across the tropics (DeMott et al., 2015;Li et al., 2013).During suppressed MJO conditions, reduced cloud cover leads to increased solar radiation at the ocean surface and decreased wind-driven mixing in the surface ocean mixed layer.The surface mixed layer, with typical depths of 10-100 m in the tropics (de Boyer Montégut et al., 2004), is characterised by nearly uniform profiles of temperature and salinity.During the reduced mixing conditions associated with the suppressed phase of the MJO, the mixed layer shoals, leading to a reduction of its heat capacity and enhancing the increase in sea-surface temperatures (SSTs) associated with the increased solar radiation.Warm SSTs moisten the low-level atmosphere via increased evaporation, creating atmospheric instability and promoting convection.During enhanced convective MJO conditions, lower incoming solar radiation at the ocean surface and increased upward latent heat flux due to strong surface winds leads to increased mechanical mixing, cooling of SST, and deepening of the ocean surface mixed layer (Drushka et al., 2012;Liu et al., 2021).
Diurnal changes in solar radiation and surface winds can lead to the development of diurnal warm layers a few metres deep, superimposed on the deeper, residual mixed layer.Seaglider observations during suppressed MJO conditions show that such layers can increase the temperature of the top few metres in the ocean by 0.8 • C, with a daily mean increase of 0.2 • C (Matthews et al., 2014).The diurnal warm layer reaches a maximum at approximately 1500 h local solar time (LST), and then disappears overnight due to nocturnal mixing.In line with observations, increased diurnal variability of SST is found to increase mean SST in the tropics in coupled model simulations (Bernie et al., 2008;Hsu et al., 2019;Seo et al., 2014).Large diurnal variability of SST can lead to increased specific humidity in the lower troposphere, affecting some simulated MJO events (Hsu et al., 2019).Increased vertical resolution in the upper ocean has also been shown to improve MJO predictions (Ge et al., 2017;Ma & Jiang, 2021).More frequent coupling in coupled models can also lead to stronger diurnal variability of SST and better onset and evolution of MJO convection (Seo et al., 2014).The cooling in the upper ocean due to the passage of the MJO in coupled simulations can also lead to improved eastward MJO propagation (Savarin & Chen, 2022).Accurate representation of two-way feedbacks between the upper ocean and atmospheric convection is essential for future improvements in MJO predictions.
Whereas coupled ocean-atmosphere climate models are widely used, there are only a few operational short-to extended-range NWP systems that use coupled configurations (Mogensen et al., 2017;Smith et al., 2018).Recently, the UK Met Office has developed a deterministic coupled ocean-atmosphere model, running in real time since 2016, alongside the atmosphere-only version of the model.Early hindcast experiments of Shelly et al. (2014) showed that the coupled model outperformed the atmosphere-only model during two strong MJO events in 2009 and 2010.However, the results of Vellinga et al. (2020) over three boreal winters showed little difference between the models in terms of MJO prediction skill.The main difference between the models was found in MJO propagation speed, with the coupled model predicting faster MJO propagation than the atmosphere-only model.In this article, we examine MJO performance using 5 years of data from real-time coupled and atmosphere-only NWP systems of the UK Met Office, expanding on the study of Vellinga et al. (2020), and using process-based diagnostics to determine the mechanisms that lead to the different MJO simulations in the coupled and atmosphere-only models.The model specifications, data, and methodology are described in Section 2. In Section 3 we present general MJO performance for both models and mechanisms leading to performance differences between the coupled and atmosphere-only models.Discussion and conclusions follow in Section 4.

Model specifications
The data used in this study were simulated with coupled ocean-atmosphere and atmosphere-only NWP systems of the UK Met Office running daily since May 1, 2016.

Global sea ice (GSI) version
May The atmosphere-only model was the operational forecast model at the time at the Met Office, and the ocean component of the coupled model was the operational ocean forecast model.Models were initialised at 0000 UTC in real time out to 10-and 7-day lead times for the coupled and the atmosphere-only models respectively.Both models yield 1,857 forecast initialisations between May 1, 2016, and May 31, 2021.The models use the same atmosphere and land components, with the addition of ocean and sea-ice models for the coupled version.Table 1 shows a summary of changes in resolution, number of vertical levels, and model components (and their references) that occurred during this study.
Both models use a mass-flux convection scheme (Gregory & Allen, 1991;Gregory & Rowntree, 1990) that allows shallow, mid-level, and deep convection.For the first 15 months of the data period, the atmosphere and land components used a horizontal resolution of N768 (0.2348 • longitude and 0.1568 • latitude).From July 12, 2017, the models were upgraded to N1280 (0.148 • longitude and 0.098 • latitude; ∼15 km and ∼10 km at the Equator).Prior to September 2018, the coupled model used an extra 15 vertical levels in the stratosphere, later changed to match the atmosphere-only model number of levels.The ocean component of the coupled model uses the Nucleus for European Modelling of the Ocean (NEMO) consortium ocean model (Madec et al., 2017), with horizontal resolution of 0.25 • and 75 vertical levels, eight of which are in the top 10 m of the ocean.The ocean-sea ice and atmosphere-land components are initialised with uncoupled data assimilation systems.The atmosphere-land component uses a four-dimensional variational data assimilation system (Rawlins et al., 2007) (hereafter, "UM analysis"), initialised at 0000 UTC, with SST and sea-ice concentrations from the Operational Sea Surface Temperature and Ice Analysis (OSTIA) (Donlon et al., 2012) assimilation system, updated by Fiedler et al. (2019) and Good et al. (2020).The initial SST and sea-ice concentrations are held constant throughout the atmosphere-only forecasts.The ocean-sea-ice component uses the Forecast Ocean Assimilating Model (FOAM)-NEMOVAR data assimilation system from Blockley et al. (2014) and Waters et al. (2015) (hereafter, "FOAM").The coupled model exchanges information between ocean-sea ice and atmosphere-land components every 1 hr.For more detailed descriptions of model configurations, see Vellinga et al. (2020, Section 2).

Real-time multivariate MJO index
MJO performance is quantified using the real-time multivariate MJO (RMM) index, originally from Wheeler and Hendon (2004).Full methodology on how the indices are calculated can be found in Gottschalck et al. (2010) and the references therein.The index uses daily anomalies of top-of-atmosphere outgoing long-wave radiation (OLR) and zonal winds at 850 and 200 hPa.The indices RMM1 and RMM2 represent the principal component time series of the dominant spatial structures (empirical orthogonal functions) of the data.The combination of RMM1 and RMM2 defines eight MJO phases depending on the location of the MJO convection in the tropics, with phases 8 and 1 being in the Western Hemisphere and Africa, phases 2 and 3 in the Indian Ocean, phases 3 and 4 in the MC, and phases 6 and 7 in the western Pacific.The amplitude of the MJO is defined as (RMM1 2 + RMM2 2 ) 1∕2 .Here, the active MJO is defined by days when amplitude ≥1.0.We use two RMM indices for model verification: the Wheeler-Hendon index (Wheeler & Hendon, 2004; retrieved from http://www.bom.gov.au/climate/mjo) and RMM indices calculated from the UM analysis, using daily means from runs initialised at 0000 UTC and 1200 UTC.RMM indices for both models are calculated from runs initialised at 0000 UTC.The model indices are compared with these two datasets using four standard scalar metrics following Lin et al. (2008) and Rashid et al. (2011): bivariate anomaly correlation coefficient, root-mean-square error (RMSE), amplitude error, and phase error.The bivariate anomaly correlation coefficient corresponds to the spatial correlation between forecasts and observations.A model is considered skilful when RMSE < √ 2 and correlation > 0.5 (Lin et al., 2008).Amplitude error is negative (positive) when the model underestimates (overestimates) MJO RMM amplitude.Phase error in the (RMM1, RMM2) plane is defined as an angle (in degrees) and is positive (negative) when the MJO is ahead of (behind) the observations.The RMM statistics are calculated for the boreal winter season (November-April) for active MJO days between May 1, 2016, and May 31, 2021.The same analysis was performed for all available data and winter season only data.Qualitatively, no notable difference was observed in RMM skill metrics between these two periods.

Composites
Composite maps are calculated from daily means of meteorological variables.The high-resolution original model data are regridded to N180 (1 Subsurface ocean data were processed along an equatorial transect to study the vertical profile of ocean-atmosphere interaction in the coupled model.The mixed layer depth is defined following Drushka et al. (2014), as the depth where the potential density  change from the potential density at a reference depth of 10 m is greater than a threshold given by where T ref and S ref are the temperature and salinity at the reference depth 10 m, P 0 is surface pressure, and ΔT = 0.8 • C is chosen as the optimal value following Kara et al. (2000).The reference depth of 10 m was chosen deliberately to remove the effects of the diurnal cycle of temperature on mixed layer depth (e.g., de Boyer Montégut et al., 2007;Hosoda et al., 2010).Temperature and salinity from the coupled model were interpolated in depth to every 1 m resolution between 0 and 1,000 m before calculating mixed-layer depth.Potential density calculations were obtained using Python package gsw v3.4.0 based on definitions from Gibbs SeaWater Oceanographic Toolbox of TEOS-10 (McDougall & Barker, 2011).The coupled model diurnal warming (dSST) is defined as the difference between 1500 h and 0600 h Local Solar Time SST.

Observational datasets
Observed daily interpolated OLR was obtained from the National Oceanic and Atmospheric Administration at 2.5 • × 2.5 • resolution (Liebmann & Smith, 1996).ERA5 reanalysis data was used for hourly 10 m wind speed (Hersbach et al., 2020) at 0.25 • × 0.25 • resolution.The National Oceanic and Atmospheric Administration OLR and the European Centre for Medium-Range Weather Forecasts reanalysis v5 (ERA5) winds were interpolated onto a 1 • × 1 • grid for comparison with model data.Daily mean short-wave (SW) radiation was obtained from the CERES SYN1deg dataset at 1 • × 1 • resolution (Rutan et al., 2015).

MJO model performance
Both the coupled and the atmosphere-only models are skilful in predicting the MJO, with bivariate correlation coefficients above 0.94 within the first seven lead days, reaching just above 0.88 for the coupled model by lead day 10 (Figure 1a).Regardless of the analysis dataset used to compare the models, there is little difference between the models in the bivariate correlation coefficients.The RMSE for Wheeler-Hendon indices (Figure 1b, dashed lines) is twice as large as the UM analysis RMSE on lead day 1 (Figure 1b, solid lines).As both models are initialised from the UM analysis, the initial RMSE is expected to be low in this case, but it converges towards the Wheeler-Hendon RMSE later in the forecast.Overall, both models are within the skilful RMSE threshold, reaching 0.62 by lead day 7.
The extended forecasts from the coupled model reach 0.90 RMSE by lead day 10, still within the threshold.
The MJO is too weak in both models at all lead times (Figure 1c), with smaller amplitude errors for the coupled model than the atmosphere-only model at lead day 3 and beyond.
The largest difference between the models is observed in phase error (Figure 1d).At lead day 1, both models predict the MJO to the east of the verification datasets.Afterwards, the coupled model tends to predict the MJO further east, with a linear increase in phase error reaching >5 • (in RMM phase space) by lead day 10.The atmosphere-only model predicts the MJO further west than the verification datasets at lead day 2 and beyond, with a constant phase error at around −1.5 • (in RMM phase space).This implies correct MJO propagation speed in the atmosphere-only model from lead day 2, albeit with the MJO anomalies placed too far to the west.Phase error in the coupled model linearly increases at a rate ∼0.6 • per day (in RMM phase space) compared with the atmosphere-only model (Figure 2a).The positive phase error difference can be interpreted as the coupled model MJO anomalies located to the east of the MJO anomalies simulated by the atmosphere-only model.During the study period, the average RMM phase speed of the MJO was 5.2 • per day (and 4.9 • per day for MJO events that stayed active crossing into the MC, i.e.RMM phase 5).The increase in phase angle error between the models of ∼0.6 • per day (in RMM phase space) is equivalent to ∼12% per day increase in MJO phase speed in the coupled model compared with the atmosphere-only model.
The same forecast skill analysis was performed for forecasts split by initial MJO phase (from phases 1 to 8 defined by the Wheeler-Hendon indices, hereafter referred to as MJO phases).The bivariate correlation coefficient (Figure 3a,b) and RMSE (Figure 3c,d) are similar between the models within seven lead days across all phases.Forecasts initialised in phases 4-6 perform better than forecasts initialised in other phases.In particular, forecasts starting in phases 1-3 show skill that drops below 0.88 after lead day 8 for the coupled model.During MJO phases 1-3, the enhanced convection is present over the Indian Ocean.By lead day 8 the convective envelope moves eastward, reaching the MC.The MC is known to produce a so-called "barrier effect" that leads to a weaker MJO or its total disappearance (e.g., Zhang & Ling, 2017).The barrier effect tends to be stronger in models than observed (Kim et al., 2014;2019;Liu et al., 2017;Seo et al., 2009;Vitart, 2017;Xiang et al., 2015).The decrease in coupled model performance at lead day 8 and beyond in initial MJO phases 1-3 is likely a result of the barrier effect present in the model.
Variations in amplitude error are larger when individual phases are considered (Figure 3e,f).Generally, amplitude is underestimated in phases 5-8 up to lead day 7 and in phases 1-3 after lead day 7.This coincides with active MJO convection over the MC, where both models display weaker convection and circulation than observed.The phase error shows that the atmosphere-only model simulates the MJO too far to the west across most phases (Figure 3h), consistent with the slow propagation seen in Figure 1d.The coupled model simulates the MJO too far to the east across all MJO phases at lead day 7 and beyond (Figure 3h).The phase error difference between the coupled and atmosphere-only models in phases 1-3 displays a similar trend to the average in Figure 2a, with ∼0.6 • increase per day in RMM phase space (Figure 2b).In phase 4, the increasing trend is observed up to lead day 4. Afterwards, the phase error difference between the models is steady and positive, implying that both models simulate MJO at a similar speed, with the coupled model anomalies located to the east of the atmosphere-only model anomalies.In initial phases 5-8, the coupled model MJO anomalies are located to the west of the atmosphere-only MJO anomalies until lead day 3, contrary to the average trend.Afterwards, the phase error difference linearly increases between the models, consistent with the average behaviour in Figure 2a.

MJO convection in the models and SST-MJO relationship
Both models predict the MJO well with slight differences in MJO convective (OLR) anomalies by lead day 7 in different initial MJO phases.In initial MJO phase 1, the observed MJO shows enhanced convection (negative OLR anomalies) over the EIO region (5 • S-5 • N, 70 • E-90 • E) and suppressed MJO convection (positive OLR anomalies) in the western MC (Figure 4a).Both models capture this convection pattern well (Figure 4c,e).By lead day 3, the coupled model exhibits stronger convection than the atmosphere-only model in the active convective EIO region (negative OLR differences of up to 5.5 W⋅m −2 by lead day 7; Figure 4g).In the convectively suppressed central MC region (120 convection in the coupled model is less suppressed than in the atmosphere-only model from lead day 3, leading to a negative 6.5 W⋅m −2 OLR difference between the models by lead day 7 (Figure 4g).In initial MJO phase 4, observations show reversed convective anomalies to those of phase 1, with enhanced MJO convection over the EIO and the MC regions (Figure 4b).As the MJO propagates eastward in initial phase 4, the EIO becomes a region of suppressed convection from lead day 4 onward.Again, both models simulate this fairly well (Figure 4d,f).However, from lead day 3, the coupled model suppresses convection faster in the EIO region than the atmosphere-only model does, reaching a positive OLR difference of 5.9 W⋅m −2 between the models by lead day 7 (Figure 4h).
These two initial MJO phases display the largest spatial differences in OLR anomalies between the coupled and atmosphere-only model at lead day 7 (Figure 5, other phases not shown).In initial phase 1, both the EIO and central MC regions display a spatially coherent difference in OLR anomaly between the models at 95% significance level (Figure 5a).In initial phase 4, the OLR anomaly difference is at the 95% significance level across almost the entire EIO region.Since both models use the same land and atmosphere components (and hence the same cumulus parametrisation scheme), the OLR anomaly differences must be driven by different ocean boundary conditions (i.e., SSTs) in the models.
Already at lead day 1 daily mean, the models exhibit differences in MJO-associated SST anomalies1 in initial phases 1 and 4 (Figure 6a,c).These differences remain fairly constant throughout the forecast, and their evolution is consistent with the differences in OLR that develop between the models by lead day 7 (Figure 5).In initial phase 1 at lead day 1, the coupled model develops warmer MJO-associated SST anomalies in the central MC region compared with the atmosphere-only model persisted SST by 0.12 • C (Figure 6a).During the suppressed convective conditions of MJO phase 1 in the central MC (Figure 4c), warm MJO-associated SST anomalies will lead to enhanced latent heat (LH) flux into the atmosphere in the region.The evaporation linked to this LH flux exchange will moisten the low level atmosphere, and in line with the moisture mode theory for eastward MJO propagation (e.g., Sobel & Maloney, 2013), this moisture anomaly will lead to more convection in the coupled model ahead of the main MJO convective envelope.Therefore, these increased warm MJO-associated SST anomalies in the central MC region at lead day 1 will lead to increased MJO propagation in the coupled model by lead day 7.
In initial phase 4, the coupled model shows colder MJO-associated SST anomalies in the EIO region than the atmosphere-only model at lead day 1 by 0.08 • C (Figure 6c).This region is in suppressed MJO conditions in phase 4 at lead day 1 (Figure 4d), and the colder MJO-associated SST anomalies will lead to inhibited convection in the coupled model there, causing a stronger suppressed MJO phase behind the main MJO convective envelope (Figure 5b).Therefore, these SST differences between the models will consistently lead to faster propagation in the coupled model, consistent with the phase speed differences we observe in Figure 2. The next question to address is how these SST differences arise.

Diurnal warming of SSTs
In this section, the role of the diurnal cycle in SST in the coupled model is examined, as a potential explanation of the SST differences between the coupled and atmosphere-only model.

Rectification of diurnal warm layer on daily mean SST
The atmosphere-only model is initialised from the previous day's OSTIA foundation SSTs, corresponding to bulk 10 m night-time ocean temperature, which excludes the effects of diurnal warming.The coupled model initial SST is the FOAM analysis 0000 UTC instantaneous ocean temperature at the top model level at 0.51 m.Hence, the coupled model initial SSTs refer to a shallower depth than the atmosphere-only initial SSTs.Additionally, the coupled model initial SSTs will have an extra component of longitudinal variation, as incoming solar radiation depends on LST at 0000 UTC.The initialisation difference in SST between the models is, however, insignificant (not shown).
The coupled model SSTs are warmer than the atmosphere-only SSTs at lead day 1 across the tropics by 0.1-0.4• C (Figure 7a).The diurnal warming of SST in the coupled model (dSST; Figure 7b), defined here as the difference between 1500 h and 0600 h LST SST, displays a similar spatial pattern to the SST difference between the models seen in Figure 7a (spatial correlation coefficient between the two patterns is 0.65 over the warm-pool region 60-180 • E).The difference between the models is around half the magnitude of the dSST because the model difference is calculated from daily mean SSTs.The strongest dSST occurs close to the Equator, peaking in the central MC and north of New Guinea at 0.5-0.6 • C. The largest positive SST difference between the models is also recorded there, at 0.2-0. of the underlying zonal gradient in the background SST in the coupled model, leading to a larger zonal gradient in the daily mean SST in the coupled model, when compared to the foundation SST used in the atmosphere-only model.This stronger SST gradient across the Indo-Pacific warm pool could improve the propagation of the MJO in the coupled model (Hu et al., 2022).The spatial pattern of the dSST across the tropics in the coupled model (Figure 7b) is broadly consistent with moored buoy array observations (Yan et al., 2021) and reanalysis data validated with surface drifters (Bellenger & Duvel, 2009).Stronger dSST over the MC can be attributed to the minimum in surface winds in that region present in the coupled model, and ERA5 reanalysis (not shown).This is consistent with glider observations that show that weaker winds and stronger SW flux into the ocean lead to stronger diurnal warming (e.g., Matthews et al., 2014).
Further evidence of the role of the diurnal warm layer can be gained by examining profiles of temperature in the upper ocean.Whereas subdaily vertical profiles of ocean temperature are not available from the coupled model output, the difference between lead day 1 daily mean temperature and the FOAM initial condition can capture the diurnal warm-layer evolution over the warm-pool region (between 60 • E and 150 • E at the Equator).The FOAM initial condition is at 0000 UTC, corresponding to 0400 h LST at 60 • E and 1000 h LST at 150 • E. Hence, over the warm-pool region, these initial conditions coincide approximately with the cool phase of the diurnal warm layer.The day 1 daily mean from the forecast model corresponds to the average of the diurnal cycle.Therefore, this daily mean in the warm-pool region will be warmer than the FOAM initial condition, and correspond approximately to the daily-mean strength of the diurnal warm layer.In the Western Hemisphere, the model is initialised during the peak diurnal warm-layer strength; therefore, the daily mean temperature profile in this region will be colder than the initial condition.
The coupled model ocean temperature shows an increase of >0.1 • C in the top 5 m of the Indo-Pacific Ocean basin from the initial condition to the day 1 daily mean (Figure 8).Crucially, the warming is not uniform across the region; the strongest and deepest warming occurs over the MC (0.4-0.6 • C), whereas in the EIO region the warming is weaker at 0.1-0.3• C.This is consistent with the patterns of SST difference between the models and the coupled model dSST in Figure 7a,b.Subsequent daily mean changes (after day 1) in warm-pool SST and upper ocean temperature in the coupled model are much smaller (less than 0.1 • C; not shown).We also note that the Western Hemisphere records a cooling in the upper 5 m of similar magnitude to the warming seen in Figure 8; local time in the Western Hemisphere at this time corresponds to the cooling phase of the diurnal warm layer.Hence, the coupled model resolves the diurnal warm-layer formation in the upper 5 m of the ocean model.This process leads to elevated daily mean SST in the coupled model and contributes to faster eastward MJO propagation via surface flux exchange.

Diurnal warm-layer strength dependence on MJO phase
The systematic warming in the coupled model SSTs at lead day 1 compared with the atmosphere-only model SSTs is present across all MJO phases.It has magnitude ∼0.2 • C in the EIO region (Figure 9a) and ∼0.3 • C in the central MC region (Figure 9b), and this magnitude is approximately constant across all MJO phases.However, there is considerable variation across MJO phases in the MJO-associated SST anomalies after the removal of the mean and annual cycle, and 20-to 200-day bandapass temporal filtering applied (Figure 9c,d).These MJO-associated SST anomalies will lead to MJO convection anomalies (Matthews, 2004;Woolnough et al., 2001).Moored buoy array observations show that diurnal warming indirectly rectifies MJO-associated SSTs (Yan et al., 2021).The atmosphere-only MJO-associated SST anomalies at lead day 1 roughly correspond to the OSTIA dataset of MJO-associated SSTs anomalies 2 and by definition exclude any diurnal warming effects.The coupled model resolves the diurnal warm-layer formation; therefore, according to findings of Yan et al. (2021)  The dSST is indicative of diurnal warm-layer strength in the ocean, as seen in Figures 7b and 8.This layer develops on days with weak surface winds and strong incoming solar radiation conditions.In the EIO region, the strongest diurnal warming is observed in phase 8 (Figure 9e) when low 10 m wind speed (Figure 10c) and high surface SW radiation flux into the ocean occur (Figure 10a).In phase 4, however, this region experiences the highest winds and moderate SW flux into the ocean, resulting in the weakest dSST of all MJO phases here.The same pattern can be observed in the central MC, where the largest dSST and largest warm MJO-associated SST anomalies occur during phase 1 (Figure 10f,d), with the lowest 10 m wind speeds (Figure 10d) and highest SW flux into the ocean (Figure 10b).
The pattern of anomalous dSST in the coupled model at lead day 1 also correlates well spatially with the difference in MJO-associated SST anomalies between the models (Figure 6a,b for initial MJO phase 1, and Figure 6c,d for initial MJO phase 4; other MJO phases not shown).The magnitude of the anomalous dSST in each MJO phase appears to explain the majority of the difference in the MJO-associated SST anomaly between the coupled and atmosphere-only model at lead day 1.At lead day 1, the atmosphere-only model MJO-associated SST anomalies can be roughly regarded as the OSTIA dataset of non-diurnally resolving MJO-associated SST anomalies.At longer lead times, both the coupled model and the OSTIA MJO-associated SST anomalies evolve, albeit that any cooling or warming recorded in the coupled model is stronger than OSTIA (not shown).The majority of the difference between the two on a subweekly time-scale can be explained by the magnitude of the anomalous dSST in the coupled model (not shown).The diurnal warm-layer formation in the coupled model is therefore the main process that modulates the MJO-associated SST anomalies at lead day 1.The MJO conditions at lead day 1 enhance or suppress the strength of the diurnal warm layer in the coupled model, rectifying daily mean SST and modulating the MJO-associated SST anomalies.These, in turn, lead to MJO convection differences between the models within the next seven forecast days.

Other potential sources for SST difference between models
Other potential mechanisms were considered for SST differences between the coupled and atmosphere-only model at lead day 1: surface SW radiation flux bias, mixed-layer depth variations, 10 m wind-speed bias, and latent heat flux bias.None of these were able to produce a significant magnitude change in SST.The details follow.

Surface SW flux bias
Surface SW flux biases could lead to a change in SST in the coupled model through direct heating of the ocean mixed layer.Although daily mean SW flux into the ocean is remarkably well reproduced in the EIO region across all MJO phases at lead day 1 compared with the CERES SYN1deg observed SW flux (Figure 10a), there is a systematic bias of approximately 10 W⋅m −2 in the central MC region (Figure 10b).An experimental case study was designed to test the magnitude of SW flux changes on SST Hence, for the 9.78 W⋅m −2 SW flux daily mean bias in the coupled model that was present for the January 17, 2017, control run, the expected SST increase due to the SW flux change would be 0.018 • C.However, the daily mean difference in SST between the coupled and atmosphere-only depth temperature changes (and therefore SST changes) being significantly smaller or larger respectively.Figure 10e,f shows mixed-layer depth calculated at the Equator for a reference depth of 10 m and a temperature change of ΔT = 0.8 • C, at lead day 1 for the EIO and central MC regions.There is no strong relationship between mixed-layer depth and MJO phase at lead day 1, with variations of mixed-layer depth of 10%-20% in both regions across different MJO phases.The mixed-layer temperature tendency due to surface heat fluxes over a 24 hr period can be estimated as follows: where Q net (W⋅m −2 ) is the net surface heat flux,  0 is seawater density (1025 kg⋅m −3 ), h (m) is mixed-layer depth, and C p is the specific heat of seawater (3850 J⋅kg −1 ⋅ • C −1 ).
Using the coupled model daily mean Q net at the Equator and mixed-layer depth at lead day 1, we find that changes in water column temperature over 24 hr in the coupled model are <0.02• C in the EIO region and <0.06 • C in the central MC region.These temperature tendencies are an order of magnitude smaller than the difference between the coupled and atmosphere-only model SST at lead day 1 at the Equator, accounting for <9% of the SST difference in the EIO region and 11%-20% in the central MC region.We also note that there are no large fluctuations in mixed-layer depth that could produce substantial changes in mixed-layer temperature through entrainment of colder waters from beneath the mixed layer.
We conclude that mixed-layer processes in the coupled model are not enough to explain the majority of the SST increase in the coupled model compared with the atmosphere-only model.The coupled model also simulates a barrier layer below the mixed layer.The simulated barrier layer is of 5 m thickness and does not vary across different initial MJO phases on lead day 1 (not shown).Observations show that barrier layers thicker than 10 m, can lead to faster SST recovery post-MJO passage due to a decreased entrainment of cold water from below the barrier layer into the mixed layer (Drushka et al., 2014;Moteki et al., 2018).This mechanism is simulated in the coupled model in the EIO region in initial phase 4 at lead day 5 and beyond.The barrier layer thickens from ∼5 m at lead day 5 to ∼10 m at lead day 10 through shoaling of the mixed layer.However, this small change in the barrier-layer thickness will have a minor, secondary effect on the SST change.

3.4.3
Surface wind speed and latent heat flux biases Latent heat flux exchange between the ocean and atmosphere acts as a cooling mechanism in the upper ocean.Latent heat flux is proportional to 10 m wind speed; therefore, we calculate the fractional wind speed bias in the coupled model compared with the ERA5 and multiply it by the coupled model surface latent heat flux to obtain the latent heat flux bias that would result from the coupled model wind speed bias.Equation (3) can be used then to calculate the latent heat flux contribution to the mixed-layer temperature tendency, with Q net replaced by latent heat flux bias.The coupled model has generally stronger winds than ERA5 by up to 7% and 15% in the EIO and central MC regions respectively at lead day 1 at the Equator (Figure 10c,d).Stronger winds will lead to an increase in latent heat flux into the atmosphere and increase the cooling effect in the upper ocean.In the EIO region, this bias yields a negative temperature tendency up to −0.005 • C for MJO phases 2-8, with MJO phase 1 yielding a positive temperature tendency of 0.001 • C. In the central MC region, the temperature tendency due to wind speed bias is negative for all MJO phases, up to −0.01 • C. Therefore, this latent heat flux bias due to the wind speed bias will lead to a slight cooling in the coupled model and cannot explain the observed increase in SST compared with the atmosphere-only model.

DISCUSSION AND CONCLUSIONS
The coupled ocean-atmosphere NWP system of the UK Met Office has been running daily since May 1, 2016.In May 2022, the coupled model was switched to an operational mode, replacing the atmosphere-only NWP system that was previously used for operational global forecasting at the Met Office.Our study reveals that the addition of an ocean model introduces new complications for the MJO forecasting in this NWP system.The inclusion of the diurnal warming of SST (dSST) in the coupled model makes it more realistic, however, its subsequent feedbacks with the MJO ultimately lead to a stronger than desired increase in the MJO propagation speed.
Both the coupled and atmosphere-only NWP models of the Met Office predict the MJO skilfully out to at least 7 lead days.However, the coupled model simulates faster eastward MJO propagation than the atmosphere-only model by 12% (in RMM phase space).This increase is caused by the coupled model's ability to resolve the diurnal warm-layer formation in the upper ocean, the effects of which are not present in the atmosphere-only model that utilises foundation (night-time) SSTs.The dSST pattern in the coupled model correlates spatially well with the SST difference between the models.Both patterns are positive across the Indo-Pacific warm pool with peak values over the MC.This uneven distribution of dSST across the Indo-Pacific warm pool leads to a stronger SST gradient in the coupled model in the region of MJO convection.The distribution of the dSST in the coupled model across the tropics is broadly consistent with observations (Bellenger & Duvel, 2009;Yan et al., 2021).Glider observations also show that weak surface winds and large SW flux lead to stronger diurnal warming (Matthews et al., 2014).Consistently, the strongest dSST in the coupled model is recorded over the MC, coinciding with the weakest surface winds in the model, and ERA5.An accurate representation of the SST gradient across the Indian Ocean is found to favour more moisture to the east of the MJO convection in coupled model simulations and leads to better eastward propagation of the MJO (Hu et al., 2022).Hence, this stronger SST gradient in the Met Office's coupled model could lead to more coherent MJO propagation.
Observations show that the dSST in the Indo-Pacific warm pool rectifies onto the intraseasonal SSTs (Yan et al., 2021).Strong dSST leads to more moist static energy ahead of the active MJO phase, leading to an earlier onset of MJO convection (Itterly et al., 2021).This feedback is simulated by the Met Office's coupled model.The MJO conditions in the coupled model dictate the strength of the dSST.The dSST rectifies onto the MJO-associated SST anomalies and those anomalies feed back into the MJO convection within the next few forecast days and lead to MJO propagation speed changes in the coupled model.This mechanism is similar to existing theories of the MJO atmosphere-ocean interaction-for a review, see DeMott et al. (2015)-and amplifies the SST patterns associated with these processes.A more detailed description of this two-way feedback between the MJO and dSST is described as follows, with a visual summary displayed in Figure 12: 1.During the convective phase of the MJO, the cloud cover is increased and wind-driven mixing is enhanced, causing the suppression of diurnal warm-layer strength (and dSST).This happens, for example, in MJO phase 4 in the EIO region.The suppressed dSST leads to colder MJO-associated SST anomalies in the coupled model at lead day 1.These colder MJO-associated SST anomalies will tend to inhibit convection in that region later in the forecast.By lead day 7, this region is in the suppressed MJO phase, to the west of the MJO convective anomalies.This mechanism will therefore result in stronger suppression behind the MJO convective envelope and lead to increased eastward MJO propagation.2. During the suppressed MJO phase, low surface winds and high incoming solar radiation lead to an enhanced dSST in the coupled model.This happens, for example, during MJO phase 1 in the central MC region located to the east of the MJO convective anomalies at lead day 1.The enhanced dSST in this region leads to warmer MJO-associated SST anomalies in the coupled model at lead day 1.These warmer MJO-associated SST anomalies will then increase moisture and convection ahead of the MJO convective anomalies, leading to faster eastward MJO propagation.
The MJO is a major source of predictability on 1to 3-week time-scales (e.g., Gottschalck et al., 2010), and therefore it is crucial for models to predict it well.Atmosphere-only models can skilfully predict the MJO out to 10 lead days (Woolnough et al., 2007), but over longer lead times the coupled NWP systems tend to outperform the atmosphere-only ones (e.g., Kim et al., 2014;Vitart, 2017).There is still much room for improvement, as models generally tend to simulate an MJO that erroneously decreases in propagation speed with lead time (Vitart, 2017).Previous studies showed that diurnal variability of SST in coupled models can lead to improved MJO predictions and simulations (Bernie et al., 2007;2008;Ge et al., 2017;Hong et al., 2017;Seo et al., 2014;Tseng et al., 2015).High vertical resolution near the ocean surface is found to be the key to stronger dSST in coupled models and increased intraseasonal SST variability that feeds into MJO convection (Ge et al., 2017;Hsu et al., 2019;Tseng et al., 2015).Consistent with these studies, we show that the dSST plays a crucial role in representing the MJO and should be considered for model improvements, especially for those models that struggle to simulate eastward propagation of the MJO.
Particular care should be taken with different SST datasets used by models.Foundation SSTs are often used for comparisons with model-simulated SSTs.Such a comparison could lead to misleading conclusions: a naive analysis could suggest that the Met Office's coupled model needs to correct the "warm bias" in the tropics.Instead, this "warm bias" is a manifestation of a real mechanism, the diurnal layer formation in the upper ocean, that is not represented in the OSTIA dataset.This mechanism may, however, be too strong in the coupled model, leading to erroneously fast MJO propagation speed.Coupled models that lack high vertical resolution near the ocean surface could potentially benefit from parametrising this mechanism.Many climate models from the Coupled Model Intercomparison Project Phase 6 (CMIP6) have the top ocean model level thickness larger than 5 m (see Wang et al., 2022, Table 1) and hence are unlikely to accurately simulate diurnal variability in the upper ocean.
The Met Office uses different configurations of the same model for weather predictions and climate simulations.The too-fast eastward-propagating MJO is also present in the seasonal configuration of this coupled NWP system (∼5 • phase error for the first 4 weeks of the forecast; Vitart, 2017) and in the CMIP6 climate configuration (HadGEM3 model; Ahn et al., 2020).Both the seasonal forecast system and the climate model are at a lower atmospheric horizontal resolution (N216) than the coupled NWP model analysed here.All three configurations, however, use the same vertical and horizontal resolutions in the ocean component of the model.It is likely that the two-way feedback between the MJO and the diurnal warm layer is present across all those configurations, irrespective of the atmospheric horizontal resolution.The high-resolution coupled NWP system of the Met Office predicts the MJO skilfully to at least 10 forecast days; however, the too-fast propagating MJO may present challenges for weather predictions past 2 weeks, and for longer term climate projections.The increase in MJO speed in the coupled model can lead to faster onset of teleconnection patterns; for example, the CMIP6 HadGEM3 model is found to underestimate the North Atlantic Oscillation response to the MJO (Skinner et al., 2022).The North Atlantic oscillation is a key component of northern Europe variability; thus, improving the MJO will improve seasonal predictions and climate projections over the UK.Lastly, the next generations of coupled models will be at higher atmospheric horizontal resolution and ultimately convection permitting.Our findings demonstrate the importance of investigating how the diurnal warm layer manifests in these models and the subsequent effects on the MJO.

F
I G U R E 1 RMM (real-time multivariate Madden-Julian Oscillation [MJO]) skill statistics as a function of lead day.(a) Bivariate correlation coefficient; (b) root-mean-square error; (c) amplitude error; (d) phase error.Daily coupled (blue) and atmosphere-only model (red) data are compared for boreal winter season (November-April) and active MJO days only with UM analysis (solid) and Wheeler-Hendon indices (dashed).

F
I G U R E 2 Phase error difference between the coupled and atmosphere-only model for (a) all Madden-Julian Oscillation (MJO) phases combined and (b) split by the initial real-time multivariate MJO phase that the forecast started in (solid contours for positive error and dashed contours for negative error).

F
I G U R E 3 RMM (real-time multivariate Madden-Julian Oscillation [MJO]) skill statistics as a function of lead day split by the initial RMM phase that the forecast started in.(a, b) Bivariate correlation coefficient; (c, d) root-mean-square error; (e, f) amplitude error; (g, h) phase error.Solid contours for positive values and dashed contours for negative values.

F
I G U R E 4 Hovmöller diagrams of daily mean composites of 20-to 200-day filtered boreal winter outgoing long-wave radiation (OLR) anomaly averaged over the equatorial band (5 • S-5 • N), for forecasts initialised in Madden-Julian Oscillation (MJO) phases 1 and 4. (a, b) Observed; (c, d) coupled model; (e, f) atmosphere-only model; (g, h) difference between coupled and atmosphere-only models.Vertical dashed lines represent equatorial Indian Ocean and central Maritime Continent regions.Initially active MJO forecasts only.Number n denotes the amount of independent events used in the composite (total number of days used displayed in the parentheses).Solid contours for positive values and dashed contours for negative values.

F
I G U R E 5 Difference at lead day 7 between composite daily means of coupled and atmosphere-only model 20-to 200-day filtered boreal winter anomaly of outgoing long-wave radiation (OLR) for forecasts initialised in Madden-Julian Oscillation (MJO) phases (a) 1 and (b) 4. Initially active MJO forecasts only.Number n denotes the amount of independent events used in the composite (total number of days used is displayed in the parentheses).Real-time multivariate Madden-Julian Oscillation (MJO) indices from Wheeler-Hendon.EIO: equatorial Indian Ocean; MC: Maritime Continent.The yellow contour outlines differences significant at the 95% level.
3 • C. The dSST is weaker in the Indian Ocean basin (up to 0.5 • C), where a weaker SST difference also occurs, albeit still positive at 0.1-0.3• C. Hence, the zonal gradient in dSST compounds the effects F I G U R E 6 Composite daily mean at lead day 1 for (a) coupled minus atmosphere-only model Madden-Julian Oscillation (MJO)-associated sea-surface temperature (SST) anomaly difference in MJO phase 1; (b) anomalous coupled model diurnal warming of SST in MJO phase 1; (c) coupled minus atmosphere-only model MJO-associated SST anomaly difference in MJO phase 4; (d) anomalous coupled model diurnal warming of SST in MJO phase 4. Boreal winter and initially active MJO forecasts only.

F
I G U R E 7 Composite daily mean at lead day 1 for (a) coupled and atmosphere-only model sea-surface temperature (SST) difference for all Madden-Julian Oscillation (MJO) phases; (b) coupled model 1500 h and 0600 h local solar time SST difference for all MJO phases.Boreal winter and initially active MJO forecasts only.Solid contours for positive values and dashed contours for negative values.
, the dSST in the coupled model may indirectly affect its MJO-associated SST anomalies.Indeed the coupled model dSST values (Figure 9e,f) correlate well with the MJO-associated SST anomalies in the 2 The atmosphere-only model in operational mode uses the previous day's OSTIA SSTs as the initial condition.Therefore, the persisted SSTs used by the atmosphere-only model are OSTIA dataset lagged by 1 day.

F
Composite mean vertical section of ocean temperature change, from Forecast Ocean Assimilating Model initial condition to lead day 1 mean (centred at 12 hr lead time), from the coupled model at the Equator for all Madden-Julian Oscillation (MJO) phases.Boreal winter and initially active MJO forecasts only.Model levels are displayed as black dots at 180 • E. Solid contours for positive values and dashed contours for negative values.coupled model (blue lines in Figure 9c,d); the coldest MJO-associated SST anomalies in the coupled model in both regions occur during the time of the weakest dSST (around MJO phases 3-5); and conversely, the warmest MJO-associated SST anomalies are present when the coupled model exhibits the strongest dSST (around MJO phases 8-2).

F I G U R E 9
Composite lead day 1 daily means for coupled (blue) and atmosphere-only (red) models for (a, b) sea-surface temperatures (SST), (c, d) Madden-Julian oscillation (MJO)-associated SST anomalies, and (e, f) diurnal warming of SST as a difference between 1500 h and 0600 h local solar time SST.Equatorial Indian Ocean (EIO) and central Maritime Continent (MC) region extents in Figure 5. Boreal winter and initially active MJO forecasts only. in the coupled model in the first 24 hr of the forecast.Forecasts were initialised on January 17, 2017, when an active MJO was in phase 1 (suppressed convection conditions in the central MC region).The experiments involved changing top-of-the-atmosphere radiation subgrid variability of cloud water content and tuning shallow cumulus clouds ("cloud erosion parameter") to artificially force a change in surface SW flux.Seven experiments showed a linear relationship between SW flux and SST daily mean changes from the control coupled run (Figure 11):

F
I G U R E 10 Composite lead day 1 daily means for coupled (blue) and atmosphere-only (red) models for: (a, b) downward short-wave (SW) surface flux, with observed values from CERES SYN1deg (Rutan et al., 2015); (c, d) 10 m wind speed, with observed values from European Centre for Medium-Range Weather Forecasts reanalysis v5 (ERA5) (Hersbach et al., 2020); (e, f) mixed-layer depth at the Equator for reference depth 10 m and ΔT = 0.8 • C. Boreal winter and initially active Madden-Julian Oscillation (MJO) forecasts only.model in the central MC region on January 17, 2017, stands at 0.31 • C (Figure 11), an order of magnitude larger than what the linear regression suggests.The coupled model response to this SW flux bias is therefore not large enough to explain the much larger difference in SST between the models.3.4.2 Mixed-layer depth variations in the coupled model A large change in mixed-layer depth in the coupled model would lead to surface fluxes being distributed over a greater or smaller depth, and subsequent mixed-layer (W•m -1 ) F I G U R E 11 Lead day 1 daily mean sea-surface temperatures (SST) sensitivity in central Maritime Continent region in a coupled forecast initialisation on January 17, 2017, to downward surface short-wave (SW) flux perturbations achieved with varying cloud erosion parameter and subgrid variability of cloud water content at the top of the atmosphere.See region extent in Figure 5.

*
Compared with a model that does not resolve the diurnal warm layer Schematic diagram of the Madden-Julian Oscillation (MJO) modulation of the diurnal warm-layer strength in the upper ocean and its subsequent rectification of daily mean sea-surface temperatures (SST) and MJO-associated SST anomalies, leading to faster eastward MJO propagation within seven forecast days.During enhanced MJO conditions, the diurnal warm layer is suppressed, leading to a colder MJO-associated SST anomaly in a model that resolves the diurnal cycle of SST.This colder anomaly will lead to decreased latent heat flux into the atmosphere and stronger suppression of MJO convection within the next few forecast days.During suppressed MJO conditions, the diurnal warm layer is enhanced, leading to a warmer MJO-associated SST anomaly in a model that resolves the diurnal cycle of SST.This warmer anomaly will lead to enhanced latent heat flux into the atmosphere, leading to more convection ahead of the MJO convective anomalies.Both mechanisms will lead to faster eastward MJO propagation in a model that resolves the diurnal warm layer compared with one that does not.Left panel modified afterYan et al. (2021).
Model specifications summary.
TA B L E 1