2.1 Model description

In this study, we have conducted four experiments, NEW-assim, NEW, OLD-assim, and OLD. In NEW-assim and NEW, we used the MIROC-ES2L (Hajima et al., 2020), and in OLD-assim and OLD, we used the MIROC-ESM (Watanabe et al., 2011). The former is newly developed for CMIP Phase 6 (CMIP6; Eyring et al., 2016), whereas the latter is an official model of CMIP5. The physical core model of MIROC-ES2L is MIROC5.2, which is a minor update of MIROC5 (Watanabe et al., 2010; Tatebe et al., 2018), whereas the physical core model of MIROC-ESM is MIROC3m (K-1 model developers, 2004). The horizontal resolution of the atmospheric component of MIROC-ES2L (MIROC-ESM) has T42 spectral truncation (i.e., approximately 300 km) with 40 (80) vertical levels up to 3 hPa (0.003 hPa). The oceanic component of MIROC-ES2L has a horizontal tripolar coordinate system. In the spherical coordinate portion south of 63^∘ N, the longitudinal grid spacing is 1^∘, while the meridional grid spacing varies from approximately 0.5^∘ near the Equator to 1^∘ in midlatitude regions. There are 62 vertical levels in a hybrid σ–z coordinate system, the lowermost of which is located at a depth of 6300 m. The oceanic component of MIROC-ESM has a horizontal bipolar coordinate system: the longitudinal grid spacing of the oceanic component is approximately 1.4^∘, while the latitudinal grid intervals vary gradually from 0.5^∘ at the Equator to 1.7^∘ near both poles. There are 44 vertical levels in a hybrid σ–z coordinate system, the lowermost of which is located at a depth of 5300 m. The resolutions in MIROC-ES2L are higher than in MIROC-ESM. In particular, 31 (21) of the 62 (44) vertical layers in MIROC-ES2L (MIROC-ESM) are within the upper 500 m of depth. The increased number of vertical layers in MIROC-ES2L has been adopted in order to better represent the equatorial thermocline.

In NEW-assim and OLD-assim, we used the ESMs that incorporated the same simple scheme for ocean data assimilation, which comprised an incremental analysis update (IAU; Bloom et al., 1996; Huang et al., 2002). This technique is relatively simple compared with more elaborate techniques such as the ensemble Kalman filter and four-dimensional variational method, but it is widely used for decadal climate predictions (e.g., Mochizuki et al., 2010; Tatebe et al., 2012). A benefit of an IAU is its relatively low computational cost, which enables decadal- to centennial-scale integration and a variety of parameter sensitivity experiments. In an IAU, during the analysis interval from t=0 to t=τ, the governing equation including a correction term for temperature and salinity (X) that is written as follows:

where adv. is the advection term; diff. is the diffusion term; F is the surface flux term; and the final term on the right-hand side is the correction term, where α is a constant, and ΔX^a is the analysis increment. We employed the values of τ=1 d and α=0.025, and the IAU was applied at depths between the sea surface and 3000 m (Tatebe et al., 2012). The analysis increment is calculated from ΔX^a=X^a(0) −X(0), where X^a(0) is the analysis and X(0) is the model first guess at t=0; this term is kept unchanged during the analysis interval from t=0 to t=τ. For X^a(0), we used observed anomalies with respect to the observed monthly mean climatology during the 1961–2000 period. For X(0), simulated anomalies in NEW-assim (OLD-assim) with respect to monthly mean climatology in NEW (OLD) were used. Such a scheme, often called “anomaly assimilation” or “anomaly initialization”, is also used in many previous studies (e.g., Smith et al., 2007; Keenlyside et al., 2008; Pohlmann et al., 2009; Li et al., 2016, 2019; Sospedra-Alfonso and Boer, 2020). The monthly objective analysis of ocean temperature and salinity (Ishii and Kimoto, 2009) is assimilated into the model as X^a, with linear interpolation to daily data. Because observed DIC concentrations are sparse in space and time, only ocean hydrographic data are used for data assimilation in the present study. Moreover, any atmospheric observations or reanalyses are not applied.

Both NEW and OLD are the exactly same as the historical simulations designated by the CMIP6 and CMIP5 protocols, respectively, with three ensemble members for each that are bifurcated from arbitrary years of the corresponding preindustrial control simulations. The ocean data assimilation experiments, NEW-assim and OLD-assim, are bifurcated from NEW and OLD at the year 1946, respectively, and they are integrated up to the year 2005. Note that the data assimilation experiments are driven with the same external forcings as in the historical simulations. In the later sections, the model results for 1961–2005 are analyzed.

2.2 Estimating pCO₂ change at the sea surface

CO2F depends on the difference in the CO₂ partial pressure between the sea and the air, i.e.,

where pCO₂ (pCO${}_{\mathrm{2}}^{\mathrm{air}})$ is the CO₂ partial pressure in the sea (air); γ is the fraction of sea ice; and K=kα is the CO₂ gas transfer coefficient, where k represents the CO₂ gas transfer velocity (Wanninkhof, 1992, 2014), and α represents the solubility of CO₂ in seawater (Weiss, 1974). The CO₂ gas transfer velocity k is a function of wind speed and the Schmidt number (Wanninkhof, 1992). This study investigated the reproducibility of the anticorrelated relationship between CO2F and SST; therefore, the direction of the flux is important. As K does not affect the direction and the flux variation due to ENSO has larger amplitude in terms of pCO₂ than pCO${}_{\mathrm{2}}^{\mathrm{air}}$ (Dong et al., 2017), the direction of the flux is governed by the variation in pCO₂. Consequently, we evaluated the pCO₂ change at the sea surface in the equatorial Pacific.

Seawater pCO₂ values depend on temperature (T), salinity (S), the DIC concentration, and the total alkalinity (Alk); therefore, the change in pCO₂ can be expanded as follows:

where C(X) = (∂pCO${}_{\mathrm{2}}/\partial X$)ΔX (X=T, S, DIC, Alk) is the pCO₂ change due to the change in X (X=T, S, DIC, Alk), and Res., which includes second-order terms (Takahashi et al., 1993), was estimated so that the left-hand and right-hand sides in Eq. (3) are equal in this study. In Sect. 3, we evaluate the CO2F and pCO₂ variations in the equatorial Pacific in NEW-assim, NEW, OLD-assim, and OLD, and we calculate each term in Eq. (3) for each experiment.

2.3 Observation and reanalysis dataset

To assess CO2F, ocean temperature, and wind speed of the model output, we used observation or reanalysis datasets. We used SOM-FFN as the CO2F dataset (Landschützer et al., 2016, 2017, 2018). It is an estimate based on the ocean surface CO₂ observation data collection, SOCATv3 (Bakker et al., 2016), and provides monthly data since 1982. It shows significant interannual variation in CO2F in some specific regions, such as the equatorial Pacific and the high latitudes of both hemispheres (Fig. S1). In Sect. 3, we focus on the CO2F in the Niño3 region (5^∘ S–5^∘ N, 150–90^∘ W) which shows notable variation in CO2F in the equatorial Pacific. This region is also the region of maximum variability for SST (Gill, 1980). The observational COBE-SST2 was used as the SST dataset (Ishii et al., 2005; Hirahara et al., 2014). The JRA-55 reanalysis (Kobayashi et al., 2015) was used as the wind speed dataset.

Figure 1(a, c) Maps showing the correlation coefficient between monthly CO2F anomalies derived from SOM-FFN and that of (a) NEW-assim and (c) OLD-assim. The analysis period is from 1982 to 2005. The solid boxes show the Niño3 region (5^∘ S–5^∘ N, 90^∘–150^∘ W) and the dotted boxes show the Niño4 region (5^∘ S–5^∘ N, 160^∘ E–150^∘ W). (b, d) Time series of the detrended NINO3-SST (blue line) and NINO3-CO2F (red line, positive upward) anomalies simulated with (b) NEW-assim and (d) OLD-assim. Values plotted are the 1-year running mean, and the shading in panels (b) and (d) shows the ensemble spread (1σ). R denotes the correlation coefficients between the detrended ensemble mean NINO3-SST and NINO3-CO2F anomalies, with a 1-year running mean filter applied.

Download

3.1 CO₂ flux and pCO₂ anomalies in the Niño3 region

Horizontal maps of the correlation coefficients between simulated and observed CO2F values are shown in Fig. 1. The model output data were the ensemble mean and were linearly interpolated into the SOM-FFN grid. Note that the data were not detrended, and a 1-year running mean filter is applied to the monthly COF2 anomalies in the 1982–2005 period before calculating the correlation coefficients in accordance with the period for which the SOM-FFN dataset is available. CO2F in NEW-assim shows a positive correlation with SOM-FFN in the equatorial Pacific region (Fig. 1a) where significant interannual variations in CO2F are found (Fig. S1). On the other hand, CO2F in OLD-assim (Fig. 1c) is negatively correlated in the equatorial Pacific. The time series in the Niño3 region of both the 1-year running mean SST (hereafter, NINO3-SST) and CO2F (hereafter, NINO3-CO2F) anomalies simulated with NEW-assim (OLD-assim) are shown in Fig. 1b (Fig. 1d). Here, the data were detrended and monthly anomalies were calculated with respect to the 1971–2000 monthly mean climatology. The correlation coefficients between NINO3-SST and NINO3-CO2F anomalies in NEW-assim, OLD-assim, and the observations are −0.50, 0.44, and −0.75, respectively (Table 1). The results in NEW-assim are consistent with the observations, whereas those in OLD-assim are not. The correlation coefficients between the NINO3-SST and NINO3-CO2F anomalies in NEW and OLD are −0.85 and −0.67, respectively (Table 1 and Fig. S2). Note that OLD could capture the observed anticorrelated relationship between the NINO3-SST and NINO3-CO2F anomalies, but OLD-assim could not reproduce this relationship.

Table 1Correlation coefficients between the detrended 1-year running mean NINO3-SST and NINO3-CO2F anomalies in NEW-assim, NEW, OLD-assim, OLD, and the observations. The correlations coefficients in NEW-assim, NEW, OLD-assim, and OLD are for the period from 1961 to 2005 (Figs. 1 and S2), and the correlations coefficient in the observations is for the period from 1982 to 2005.

Download Print Version | Download XLSX

As the vertical direction of CO2F is determined mainly by pCO₂ at the sea surface (see Eq. 2), we further estimated each term in Eq. (3) for each model output (Fig. 2). ∂pCO${}_{\mathrm{2}}/\partial X$ in the C(X) term in Eq. (3) (X=T, S, DIC, or Alk) was estimated based on the climatological annual mean T, S, DIC, and Alk at the sea surface within the Niño3 region in each experiment. ΔX in C(X) (ΔpCO₂ on the left-hand side of Eq. 3) is the variation in X (pCO₂) associated with ENSO and was calculated by averaging the monthly mean X (pCO₂) anomalies regressed on the NINO3-SST anomalies over the entire Niño3 region. Note that the NINO3-SST anomalies are standardized by the standard deviation. In the following, we describe the anomalies during El Niño periods, whereas the opposite applies during La Niña periods. In NEW-assim, NEW, and OLD, pCO₂ decreases because the effect of the decrease in pCO₂ with decreasing DIC concentrations is larger than that of the increase in pCO₂ with warming (Fig. 2). In OLD-assim, however, the effect of the increase in pCO₂ with warming is larger than that of OLD, and the decrease in pCO₂ with decreasing DIC concentrations is smaller than that of OLD, resulting in an increase in pCO₂. As noted in Sect. 1, previous studies (Le Borgne et al., 2002; Feely et al., 2004; Doney et al., 2009a, b) have shown that the variability in the upwelling during ENSO events dominates the equatorial Pacific CO2F variations through its regulation of DIC. In the following, we discuss the temperature and vertical velocity changes associated with ENSO along the Equator.

Figure 2Changes in pCO₂ regressed onto the standardized NINO3-SST anomalies (ΔpCO₂) (dots) and its decomposition with changes in X (X=T, S, DIC, Alk), C(X), and Res. (Eq. 3) evaluated in NEW-assim, NEW, OLD-assim, and OLD. See the text for details regarding the calculation of ΔpCO₂, C(X), and Res.

Download

3.2 DIC and vertical velocity changes

A cross section of the monthly ocean temperature anomalies regressed onto the standardized monthly mean NINO3-SST anomalies along the equatorial Pacific is presented in Figs. 3 and S3 in addition to the climatological annual mean depths of the 18, 20, and 22 ^∘C isotherms. Here, monthly temperature anomalies were calculated with respect to the 1971–2000 monthly mean climatology. The observational temperature anomalies and the climatological isotherms are derived from the monthly objective analysis of ocean temperature (Ishii and Kimoto, 2009). Amplitudes of the positive (negative) equatorial temperature anomalies in the upper (lower) layer of the eastern (western) equatorial Pacific in NEW are larger than in OLD and are closer to the observations. The intensity of ENSO, defined as the standard deviation of detrended 1-year running mean NINO3-SST anomalies from 1961 to 2005, is shown in Table 2. The intensity of ENSO in NEW is estimated to be 1.17 ^∘C (Table 2), which is a bit stronger than the observed value (0.80 ^∘C). On the other hand, the intensity of ENSO in OLD is 0.43 ^∘C, which is about half as large as the observed value. In addition, the climatological mean thermocline in NEW is tighter than in OLD and is closer to the observations. The improvement in ENSO reproducibility in NEW is mainly attributed to two updates in the model configuration. The first is the implementation of an updated plume model for cumulus convection with multiple cloud types where the lateral entrainment rate varies vertically depending on the surrounding environment (Chikira and Sugiyama, 2010). The state-dependent lateral entrainment affects the strength of the convectively induced coupled air–sea processes in the eastern tropical Pacific and, thus, the ENSO amplitude in the model. More details are described in Watanabe et al. (2010). The second is the reduction of numerical diffusion due to the introduction of a highly accurate tracer advection scheme in the ocean and an increase in the vertical resolution (Prather, 1986). The equatorial thermocline in the climatic-mean state of the tropical Pacific is more diffuse in OLD than in the observations, which is partly due to numerical diffusion, especially in vertical advection (Tatebe and Hasumi, 2010), and this model bias is much alleviated in NEW. Correspondingly, the so-called “thermocline mode” (e.g., Imada and Kimoto, 2006) becomes more effective and the ENSO amplitude becomes larger in NEW. As the ENSO amplitude in NEW is larger than in OLD, the variation in the equatorial trade winds, which causes anomalous equatorial vertical velocity, is also larger in NEW.

Table 2The intensity and period of ENSO in NEW, OLD, and the observations calculated from the 1-year running mean NINO3-SST anomalies for the period from 1961 to 2005.

Download Print Version | Download XLSX

To assess the variations in zonal wind associated with ENSO, we estimated the 10 m zonal wind anomalies over the NINO4 region (5^∘ S–5^∘ N, 160^∘ E–150^∘ W; the dotted boxes in Fig. 1a) which are regressed onto the NINO3-SST anomalies (Table 3). The Niño4 region is the region of maximum variability for the equatorial trade winds (Fig. S4). Hereafter, the abovementioned regression coefficient is referred to as wind feedback (Guilyardi et al., 2006). The positive value of the wind feedback in NEW (0.92 m s⁻¹ K⁻¹) indicates westerly wind anomalies during El Niño, and this is consistent with the observational dataset, i.e., 1.02 m s⁻¹ K⁻¹. The wind feedback in OLD (0.46 m s⁻¹ K⁻¹) is about half of that in NEW and the observations, respectively.

Table 3The wind feedback and the vertical velocity feedback in NEW-assim, NEW, OLD, and OLD-assim. The wind feedback is computed as the monthly 10 m zonal wind anomalies in the Niño4 region regressed onto the monthly NINO3-SST anomalies, and the vertical velocity feedback is the monthly vertical velocity anomalies at the depth of the 20 ^∘C isotherm in the Niño3 region regressed onto the monthly NINO3-SST anomalies. The wind feedback is also evaluated from the observation dataset.

NA: not available.

Download Print Version | Download XLSX

Figure 3Anomalies of equatorial ocean temperature regressed onto the standardized NINO3-SST anomalies for NEW (a), OLD (b), and the observations (c). The contour interval is 0.1 ^∘C. The thick solid lines indicate the climatological-mean of the 18, 20, and 22 ^∘C isotherms.

Download

Cross sections of the monthly upward water velocity and DIC concentration anomalies along the Equator regressed onto the standardized NINO3-SST anomalies in NEW (OLD) are shown in Fig. 4a and c (Fig. 4b and d), respectively. By reproducing wind feedback that is consistent with the observations, the westerly wind anomalies during El Niño periods in NEW (Fig. S4c) are comparable to that of the JRA-55 reanalysis (Fig. S4i), leading to a weakening of the upward vertical velocity of approximately $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ (Fig. 4a). This weakening of the upward vertical velocity causes a decrease in the surface DIC in the eastern equatorial Pacific during El Niño periods (Fig. 4c). In OLD, the smaller wind feedback and the associated smaller westerly wind anomalies than in the JRA-55 reanalysis (Fig. S4g) lead to a weakening of the upward vertical velocity of just 10⁻⁶ m s⁻¹ in the equatorial Pacific (Fig. 4b). Although the ENSO signal in OLD is weaker than in the observations, due to a decrease in the upward vertical velocity compared with normal years, the surface DIC concentration decreases during El Niño periods (Fig. 4d). This is consistent with Dong et al. (2016) and shows that OLD is able to qualitatively reproduce the negative correlation between the SST and DIC concentration anomalies in the eastern equatorial Pacific (Fig. S2b).

Figure 4Anomalies of the equatorial vertical velocity (a, b) and DIC (c, d) regressed onto the standardized NINO3-SST anomalies for NEW (a, c) and OLD (b, d). Contour intervals are 0.5 × 10⁻⁶ m s⁻¹ in panels (a) and (b) and 2 μmol L⁻¹ in panels (c) and (d). The thick solid lines indicate the climatological mean of the 18, 20, and 22 ^∘C isotherms.

Download

Figure 5Equatorial temperature analysis increments (a, b), vertical velocity anomalies (c, d), and DIC anomalies (e, f) regressed onto the standardized NINO3-SST anomalies for NEW-assim (a, c, e) and OLD-assim (b, d, f). The contour intervals are 0.02 × 10^−6 ∘C s⁻¹ in panels (a) and (b), 0.5 × 10⁻⁶ m s⁻¹ in panels (c) and (d), and 2 µmol L⁻¹ in panels (e) and (f). The thick solid lines indicate the climatological mean of the 18, 20, and 22 ^∘C isotherms.

Download

Next, we examined the correction term in temperature due to the data assimilation, i.e., the temperature analysis increment; the final term on the right-hand side of Eq. (1); and the variations in vertical velocity and DIC concentration. Anomalies of the monthly mean temperature analysis increments, the vertical velocity, and the DIC concentration along the Equator regressed onto the standardized NINO3-SST anomalies are shown in Fig. 5. The maximum absolute value of the equatorial temperature analysis increment in NEW-assim is found at depths of 10–40 m in the eastern equatorial Pacific, which is shallower than the depth of the thermocline (Fig. 5a). In NEW-assim, the wind feedback is 0.92 m s⁻¹ K⁻¹ (Table 3), which is of the same magnitude as that in NEW (0.92 m s⁻¹ K⁻¹), and the surface wind anomalies still show a similar pattern to that in NEW (Fig. S4a–d). The westerly wind anomalies in NEW-assim lead to a weakening of the upward vertical velocity along the Equator during El Niño periods (Fig. 5c). To assess the variation in the equatorial vertical velocity associated with ENSO, we estimated the anomalies of the vertical velocity at the depth of the 20 ^∘C isotherm (the depth of the thermocline) in the Niño3 region which are then regressed onto the NINO3-SST anomalies. Hereafter, the regression coefficient is referred to as the vertical velocity feedback. The vertical velocity feedback in NEW-assim is estimated to be $-\mathrm{4.5}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹ K⁻¹, which is not significantly different from that in NEW (−3.9 × 10⁻⁷ m s⁻¹ K⁻¹) (Table 3). The negative value of the vertical velocity feedback in NEW-assim indicates a weakening of the upward vertical velocity at the depth of the thermocline during El Niño periods in the eastern equatorial Pacific (Fig. 5c). The weakening of the upward vertical velocity causes a decrease in the supply of DIC-rich subsurface water to the surface layer, leading to a reduction in the surface DIC concentration (Fig. 5e). In OLD, the temperature variations associated with ENSO at the depth of the thermocline in the eastern equatorial Pacific are smaller than observed (see Fig. 3b and c), so that the correction term forces a 0.16 × 10^−6 ∘C s⁻¹ increase in the equatorial water temperature during El Niño periods in order to realize the observed temperature variations (Figs. 5b, S3b). The wind feedback in OLD-assim is 0.48 m s⁻¹ K⁻¹ (Table 3), which is the same as in OLD, and the map of the wind speed anomalies shows a similar pattern to that of OLD (Fig. S4e–h); however, the warming due to the data assimilation procedure during El Niño periods reduces the density, leading to low-pressure anomalies. This results in anomalous cyclonic circulation and convergence and, thus, an enhancement of the upward vertical velocity at the depth of the thermocline (Fig. 5d). The vertical velocity feedback in OLD-assim is 4.1 × 10⁻⁷ m s⁻¹ K⁻¹, which has an opposite sign to that in OLD, −4.9 × 10⁻⁷ m s⁻¹ K⁻¹ (Table 3). A positive value of the vertical velocity feedback indicates the enhancement of the upward vertical velocity at the depth of the thermocline during El Niño periods, which is inconsistent with the observations. This spurious enhancement of the upward vertical velocity during El Niño periods causes an increase in the surface DIC concentration (Fig. 5f), leading to a positive correlation between the SST and CO2F (Fig. 1d), contrary to observations. We have to note here that the vertical velocity distribution (Fig. 5c) is even different from NEW in NEW-assim due to the temperature analysis increment (Fig. 4a). As already discussed, the intensity of ENSO in NEW is slightly stronger than observed (Table 2). In addition, the period of ENSO, which is defined as the peak of the power spectrum of the 1-year running mean NINO3-SST, is 5.0 years in NEW, which is longer than the 3.5 years in the observations (see Table 2). Because the ENSO characteristics in NEW are not perfectly consistent with observations, the model nature, namely the responses of the vertical velocity and DIC concentration in ENSO, are still distorted by the temperature analysis increment, even in NEW-assim. This indicates that further model improvements are needed.