National Marine Environmental Forecasting Center, Ministry of Natural Resources, Beijing 100081, China
2.
National Supercomputing Center in Wuxi, Wuxi 214315, China
3.
Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters, Nanjing University of Information Science & Technology, Nanjing 210044, China
4.
Nansen-Zhu International Research Center, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China
Funds:
The National Natural Science Foundation of China under contract Nos 42192561 and 41605052; the National Key Research and Development Program of China under contract No. 2020YFA0608804.
Subduction process is a dynamical bridge for the exchanges of heat between the atmosphere and subsurface ocean water, which is regarded as a central proxy for the ocean climate studies. Given its key indicator in climate signals, it is of importance to examine the ability of a model to simulate the global subduction rate before investigating the climate dynamics. In this paper, we evaluated the ability of 21 climate models from Coupled Model Intercomparison Project Phase 6 (CMIP6) in simulating the subduction rate. In general, the simulation ability of the models to the subduction climatology is better than that to the long-term variation trend. Based on the comprehensive analysis of climatology distribution and long-term trend of the subduction rate, GISS-E2-1-G performs better in reproducing the subduction rate climatology and IPSL-CM6A-LR can simulate positive long-term trend for both the global mean subduction rate and the lateral induction term in the Antarctic Circumpolar Current (ACC) region. However, it is still challenging to capture both the distribution characteristics of the subduction climatology and the long-term temporal trend for the 21 CMIP6 models. In addition, the model results demonstrate that, the ACC area is the major region contributing to the long-term trend of the global mean subduction rate. The analysis in this paper indicates that the poor simulation ability of reproducing the long-term trend of global mean subduction rate might be attributed to the ocean dynamics, for example, the zonal velocity at the bottom mixed layer and zonal gradient of mixed layer depth.
Subduction refers to the irreversible process in which water bodies within the ocean mixed layer penetrate downward into the main thermocline layer. This process occurs mainly during the deep winter and is an important mechanism by which the abnormal signals of the atmosphere and the ocean mixed layer can propagate to greater ocean depths. Moreover, it affects the formation of mode water to a large extent (Oka and Suga, 2005). Mode water is an important product of the subduction process. It is a water body with relatively uniform salinity, temperature and potential vorticity in the vertical direction. Mode water usually exits the upper layer of the main thermocline, and its homogeneity in physical properties contrasts sharply with the water stratification in the thermocline (Oka and Qiu, 2012). It brings the low potential vorticity water body from the ocean surface into the subsurface layer and directly affects the variation in ocean circulation and upper ocean stratification. Therefore, it has far-reaching dynamic and thermal effects on the ocean (Liu et al., 2011). In addition, mode water can reflect changes in the ocean and atmosphere during its formation. Variations in the properties, distribution, and transport processes of mode water can reflect changes in sea-air interactions, ocean heat transport, vortex activities, and circulation intensity in its cradle (e.g., Suga and Hanawa, 1995; Kubokawa, 1999; Gu and Philander, 1997; Liu and Huang, 1998; McPhaden and Zhang, 2002). In the local climate system, mode water plays a key role in the transportation of nutrients, heat, fresh water and carbon (Bates et al., 2002; Palter et al., 2005). Considering the significance of mode water in the ocean and climate and its close relationship with subduction during its formation, it has important implications for studying subduction processes.
Mode water in the world’s oceans can be divided into subtropical mode water and subpolar mode water. Thus far, mode water has been found in every basin. Subtropical mode water is closely associated with the western boundary current separation region of each subtropical circulation, while mode water has also been found in the northeast of some subtropical circulations. Subpolar mode water refers mainly to the subpolar region of the North Atlantic and the mode water to the north of the subpolar front of the Southern Hemisphere. Much work related to mode water has been done based on observations and model results. For example, many studies have demonstrated that the subduction rate of mode water is characterized by interannual and interdecadal changes (Qu and Chen, 2009; Luo et al., 2009; Liu and Huang, 2012; Ma and Lan, 2017). Other studies have examined the correlation between the production of mode water and the Pacific Decadal Oscillation (PDO) pattern (Qu and Chen, 2009; Ladd and Thompson, 2002). Subsequent studies have focused on the impact of global warming on mode water production (Luo et al., 2009, 2011). Gao et al. (2018) noted that the climate response to radiation forcing is delayed by the huge thermal capacity of the ocean under global warming. In this process, mode water, especially Subantarctic Mode Water (SAMW), plays an important role.
Many studies have reported the spatial features of subduction with observations and model simulations. Huang and Qiu (1994) estimated the subduction rate using Levitus climatology data (Levitus, 1982) and determined the spatial distribution and specific values of subduction. Subsequent studies have used various sources of data to investigate subduction in different areas (Qu et al., 2002; Suga et al., 2008; Qu and Chen, 2009; Chen et al., 2010; Oka and Qiu, 2012; Toyama et al., 2015; Ma and Lan, 2017). In addition, some researchers have examined the time variability of subduction and its possible influence. For instance, Liu and Huang (2012) analyzed the temporal variability of the global subduction rate during 1959−2006 on both interannual and decadal time scales and reported that the variation in the mixed layer depth (MLD) caused by lateral induction played a leading role in this variability. Chen et al. (2018) linked the global warming hiatus to subduction and explored the relationship between the variation in subduction and the global mean sea surface temperature; they argued that the subduction variability within the Antarctic Circumpolar Current (ACC) region has an important impact on global warming and is mainly caused by sea surface wind stress. Xia et al. (2021) used Coupled Model Intercomparison Project Phase 5 (CMIP5) projection data to explore the fast and slow responses of the SAMW in the South Indian Ocean to global warming and pointed out that the response of the SAMW is dependent both on the evolution of surface buoyancy forcing and overlying westerly wind stress.
The development of climate models, especially the implementation of the CMIP series of models, provides a good platform for us to study and compare the subduction rates in state-of-the-art climate model simulations. Subduction directly relates to the oceans and atmosphere through surface forcing and mixed layer changes. Therefore, analyzing the subduction results of CMIP6 models can give us a better understanding of the air-sea interaction simulation ability of these models, thus providing a meaningful reference for improving the corresponding physical processes in the model. The CMIP6 model outputs (Eyring et al., 2016) were recently released, but subduction simulations on a global scale in CMIP6 models have not yet been the focus of research. Therefore, a natural question arises as to how well the subduction rates in the CMIP6 models are reproduced. This study aims to use the outputs of CMIP6 historical simulations to evaluate the subduction simulations in the latest versions of climate models. Section 2 gives brief descriptions of the models, validation datasets and calculation methods. The simulation results are compared between the CMIP6 models and validation datasets in Sections 3 and 4. Section 5 provides the discussion and conclusion.
2.
Data and methods
2.1
Data
The 21 CMIP6 models used in this study are summarized in Table 1 that are available at the time of initializing this study. The monthly potential temperature, salinity and horizontal and vertical velocities of the ocean component in the historical simulations are employed. Each historical simulation is integrated from a preindustrial spin-up experiment and is forced by solar irradiation, volcanic activity, aerosol loadings, and greenhouse gas emissions for the period of 1850−2014. In this study, we select the period from 1901 to 2014 for analysis.
Table
1.
Details of the Coupled Model Intercomparison Project Phase 6 (CMIP6) models
The subduction rate climatology simulation results are compared with the gridded Argo data provided by the Asia-Pacific Data Research Center and the Simple Ocean Data Assimilation 2.2.4 (SODA 2.2.4) reanalysis (Giese and Ray, 2011). The Argo data have a horizontal resolution of 1°×1° and 27 vertical levels. The period from 2005 to 2019 is employed to calculate the subduction rate climatology. SODA monthly data including potential temperature, salinity, horizontal and vertical velocities from 1901 to 2009 are selected. In addition, the monthly surface wind from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) is used for the calculation of Ekman pumping in Argo result. The horizontal resolution is 2.5°×2.5° and the product from 2005 to 2019 is used in this study.
SODA 2.2.4 is used continuously in the trend analysis. It is chosen here due to its long assimilation time covering more than 100 a, which is closer to the time length of the CMIP6 historical simulations; thus, it provides an advantage in trend analysis in the subduction rate. For mutual conformation, the trend of Ocean Reanalysis System 3 (ORA-S3) data with a period from 1959 to 2010 is also analyzed. The model simulation and reference data are all interpolated to a 1°×1° grid.
2.2
Methods
The annual subduction rate is calculated following Huang and Qiu (1994), who tracked water particles released from the base of the ocean mixed layer in late winter for a whole year until the late winter of the following year. This method estimates the subduction rate in Lagrangian coordinates and the specific equation is
where Sann is the annual mean subduction rate; $ {w}_{\mathrm{E}\mathrm{K}} $ is the Ekman pumping velocity; hm is the mixed layer depth; $ v $ is the meridional current velocity; t1 and t2 both of which are in late winter, denote the start and end times of tracking, respectively. Late winter is defined as March in the Northern Hemisphere and September in the Southern Hemisphere. The integrated time spans a whole year T. The first term on the right-hand side of Eq. (1) represents the contribution of vertical pumping at the bottom mixed layer (BML), while the second term is the contribution from the MLD advection, which is closely associated with the horizontal distribution of the MLD and the horizontal velocity. In the CMIP6, SODA and ORAS-3 results, the vertical current velocity is conveniently obtained, therefore the interpolated vertical velocity at the base of the mixed layer is used to replace the vertical pumping term in Eq. (1), just like the calculation performed in Liu and Wu (2012) and Chen et al. (2018). The annual mean subduction represents the water mass that migrates from the surface mixed layer into the permanent pycnocline and does not return in the late winter of the following year; therefore, the annual mean subduction rate should not be negative, and negative values indicate that there is no effective detrainment, in which case we set the rate to zero. We define the MLD as the depth where the potential density is 0.125 kg/m3 larger than that at the sea surface in the vertical direction.
3.
Climatology
3.1
The mixed layer depth
The simulation skills for the MLD have a great impact on the performance of subduction rate in the models. Therefore, the reproduction of the intermittent nature of the MLD in different seasons for the CMIP6 models is checked here. In the Northern Hemisphere, spring is from April to June, summer is from July to September, autumn is from October to December, and winter is from January to March. While in the Southern Hemisphere, the seasons are opposite.
The multi-model ensemble mean (MME) of these CMIP6 models are capable of capturing the MLD spatial characteristics as revealed by Argo observations in the four seasons. The bias between the CMIP6 MME simulation and the Argo data is mainly distributed in the relatively deep MLD regions, such as the northern North Atlantic and the ACC region (Fig. 1). The MLD in these two regions are shallower than that in Argo regardless of seasons. In addition, the bias of the SODA data are quite similar.
Figure
1.
The mixed layer depth for Argo in spring (a), summer (b), autumn (c) and winter (d); e−h are the same as a−d, but for Simple Ocean Data Assimilation results; i−l are the results for multi-model ensemble mean.
Figure 2 shows the MLD simulation skills of these CMIP6 models with reference to Argo and SODA. The groups of spring and winter simulations are overall better than the other seasons according to the correlation coefficients with Argo. While the skill of summer lags behind due to the low correlation scores of CAMS-CSM1-0, FGOALS-f3-L, FGOALS-g3 and MCM-UA-1-0. For standardized deviations, most models are in the range of 0.75−1.25. In all seasons, MCM-UA-1-0 and NorESM2-LM are the most two common outliers (Fig. 2a). The result using SODA as reference data in Fig. 2b shows some differences. The amplitudes are generally large as the standardized deviations of SODA is relatively small than that of Argo. Moreover, another obvious difference is the artificial high correlation in the summer season. This may be attributed to the consistent bias for tropical MLD in both SODA and CMIP6 models, which can be seen in Fig. 3b.
Figure
2.
Taylor diagram for mixed layer depth in four seasons for the 21 Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME) with Argo (a) and Simple Ocean Data Assimilation (SODA) (b) data.
Figure
3.
Zonal averaged mixed layer depth (MLD) during the spring (a), summer (b), autumn (c) and winter (d) for the Coupled Model Intercomparison Project Phase 6 (CMIP6) models, Simple Ocean Data Assimilation (SODA) data, multi-model ensemble mean (MME), as well as the reference datasets.
Zonal averaged MLD during four seasons is analyzed in Fig. 3. The MLD shows obvious meridional and seasonal variations. There are two peaks for MLD in the latitude bands of 50°−60°S and 50°−60°N. The MLD has the maximum value in winter, followed by spring, autumn and summer. The CMIP6 models are capable in simulating the variation features for MLD. However, the biases are evidenced in the tropics in summer, and in the high latitudes band mentioned above in all seasons. For the individual models, the MLDs simulated by NorESM2-LM and CanESM5 are deeper than that derived from Argo, especially in the ACC region. MCM-UA-1-0 reproduces a much shallower MLD, also with a smaller amplitude in all latitudes. GISS-E2-1-G and CESM2-WACCM are the two outstanding pioneers in simulating the MLD in ACC region. Furthermore, GISS-E2-1-G also performs well in the middle and low latitudes which excels CESM2-WACCM in these regions. While the MLD in SODA is also shallower in the high latitude bands, the simulation in other latitudes is more consistent. It is worth mentioning that there are artificial double peaks in the summer MLD near the equator for both SODA and CMIP6 models, which is also responsible for the unreal high correlation coefficients between SODA and CMIP6 in summer season as shown in Fig. 2b.
3.2
Subduction rate
Evaluating the subduction rate climatology is one of the fundamental steps in examining the simulation ability of various models. The global subduction climatology of Argo and SODA outputs are shown in Fig. 4. The results of these two datasets share many common features in the spatial distribution, mainly in the following aspects. High subduction rate centers closely relate to the corresponding mode water regions in each ocean basin. In the Northern Hemisphere, significant subduction centers are located near the Kuroshio Extension in the western North Pacific, the eastern North Pacific and the northern North Atlantic. In the Southern Hemisphere, large subduction rates are distributed mainly at the same latitudes (30°−60°S) of the South Indian Ocean, the South Pacific Ocean and the South Atlantic Ocean, which are associated with the ACC region. The differences between Argo and SODA outputs are mainly reflected in the intensity of the subduction rate centers. In SODA result (Fig. 4b), the subduction rate is stronger in the Northwest Pacific and North Atlantic than that in Argo (Fig. 4a). However, whether in the subtropical sea area or the ACC region, the subduction derived from SODA is weaker. The three typical cases are the Northwest Pacific, the North Atlantic and ACC region, the subduction rates responding to the three regions in Argo are about 100 m/a, 400 m/a, and 400 m/a, while in SODA data, they are 200 m/a, 800 m/a, and 200 m/a, respectively. Therefore, although SODA is an assimilation data, it is competent as a reference for the CMIP6 model evaluation when the climatology is assessed.
Figure
4.
Subduction rate climatology of Argo data (a) and Simple Ocean Data Assimilation (SODA) data (b).
Figure 5 plots the global annual subduction rates for each of the CMIP6 models averaged from 1901 to 2014. The analyzed CMIP6 models are capable in simulating the subduction rate distributions as revealed by the reference data. The biases between the CMIP6 results and the observations are mainly reflected in the positions and intensity of the subduction rate centers, especially in the North Atlantic and ACC region. The subduction rate in the North Atlantic is greater in the CAMS-CSM1-0 and MPI-ESM1-2-LR, whereas weaker in CAS-ESM2-0, CESM2, CESM2-FV2, CESM2-WACCM, E3SM-1-0, and E3SM-1-1, in comparision with Argo result. Most CMIP6 models have produced narrower band for subduction rate in the ACC region, also with weaker peak values. The CAMS-CSM1-0, MPI-ESM1-2-LR and NESM3 simulate a wider and stronger band of subduction rate in the ACC region, which is closer to the observation.
Figure
5.
Subduction rate climatology for each Coupled Model Intercomparison Project Phase 6 (CMIP6) model.
The quantitative understanding of the individual model skill in capturing the subduction rate climatology is analyzed by the Taylor diagram in Fig. 6. With reference to Argo results, most models get a correlation coefficient about 0.2, while MME gets the highest score (above 0.30). For the standardized deviation, E3SM-1-0 and E3SM-1-1 get the closest deviations to Argo data; only there models reproduce a weaker amplitude, such as CanESM5, FGOALS-f3-L and GISS-E2-1-G; the remaining models all simulate a stronger amplitude of subduction rate climatology, such as CESM2-FV2 and CIESM, which are the top two obvious outliers (above 2.5).
Figure
6.
Taylor diagram of subduction rate climatology simulation of the 21 Coupled Model Intercomparison Project Phase 6 (CMIP6) models and the multi-model ensemble mean (MME) with Argo (a) and Simple Ocean Data Assimilation (SODA) (b) data.
The results with SODA data as reference (Fig. 6b) are slightly different from those of Argo. The correlation coefficients between CMIP6 models and SODA are higher, with more models around 0.3. And the number of models distributed in the high standard deviation range is also more. However, in general, the two results are quite similar, both with MME getting the highest score of correlation coefficient. In terms of standard deviations, the amplitudes of E3SM-1-0, E3SM-1-1, CanESM5, FGOALS-f3-L and GISS-E2-1-G are relatively close to those of the reference data, while the most two significant outliers are all CESM2-FV2 and CIESM. Moreover, considering the correlation coefficient and standardized deviation, GISS-E2-1-G performs better no matter with reference to Argo or SODA.
Figure 7 shows the subduction rate climatology of the CMIP6 MME and the difference relative to Argo and SODA data. In general, the MME results reproduce the distribution characteristics of the subduction rate well, and the high subduction rate centers mentioned above in the North Pacific, North Atlantic and ACC regions are all captured. Compared with Argo results, the positions of subduction rate centers shift equatorward both in western North Pacific and the ACC region. And the bias relative to SODA also suffers from this problem. The equatorward bias of subduction hot spots is consistent with the equatorward bias of the MLD in winter shown in Fig. 1d and Fig. 1l. A recent study (Hong et al., 2021) also pointed out the equatorward bias of subduction rate center and the corresponding MLD in the ACC region that may be a common problem for some CMIP6 models in simulating subduction rate hot spots.
Figure
7.
Subduction rate climatology of the Coupled Model Intercomparison Project Phase 6 (CMIP6) multi-model ensemble mean (MME) (a), the bias relative to Argo (b) and the bias relative to Simple Ocean Data Assimilation (SODA) (c).
According to Eq. (1), the subduction rate consists of two terms. One is the contribution from vertical pumping, and the other is the contribution from lateral induction associated with the slope of the mixed layer base. Figure 8 shows the climatology for vertical pumping and lateral induction derived from Argo and SODA data. In the core regions of subduction rate, the lateral induction values are obvious larger than vertical pumping and the latter mainly occurs in the subtropics. The same conclusion can be reached both in the results derived from Argo and SODA data. And the ratios between lateral induction and vertical pumping in Argo and SODA data are 1.29 and 1.37, which are quite close to each other.
Figure
8.
Vertical pumping (a) and lateral induction (b) climatology of Argo; c and d are the corresponding results of Simple Ocean Data Assimilation (SODA) data. The value in the top right-hand corner represents its contribution to subduction rate.
Figures 9 and 10 are the simulated vertical pumping and lateral induction distributions in the CMIP6 models. Some CMIP6 models are capable of capturing the relatively large contribution of lateral induction as revealed by the reference data, such as CAMS-CSM1-0, CanESM5, CAS-ESM2-0, CESM2-WACCM, CESM2-WACCM-FV2, EC-Earth3-Veg, FGOALS-f3-L, FGOALS-g3, GISS-E2-1-G, IPSL-CM6A-LR, MCM-UA-1-0, MPI-ESM1-2-LR, NESM3 and NorESM2-LM. While the contributions of vertical pumping in CESM2, CESM2-FV2, CIESM, E3SM-1-0, E3SM-1-1, FIO-ESM-2-0 and SAM0-UNICON are even greater than lateral induction which is inconsistent with that derived from reference data. The ratio of lateral induction and vertical pumping is 1.29 in Argo and the three CMIP6 models with the closest proportions are GISS-E2-1-G (1.28), NorESM2-LM (1.34), and IPSL-CM6A-LR (1.41).
Figure
9.
Vertical pumping climatology for each Coupled Model Intercomparison Project Phase 6 (CMIP6) model. The value in the top right-hand corner represents its corresponding contribution to subduction rate.
Figure
10.
Lateral induction climatology for each Coupled Model Intercomparison Project Phase 6 (CMIP6) model. The value in the top right-hand corner represents its corresponding contribution to subduction rate.
Taylor diagram for vertical pumping and lateral induction for the CMIP6 models with reference to Argo and SODA data are shown in Fig. 11. In general, the simulation skill for lateral induction is better than that for vertical pumping, which could be reflected in higher correlation coefficients. Compared to Argo results (Fig. 11a), the models’ correlation in vertical pumping simulation are almost below 0.2, only FGOALS-f3-L, FGOALS-g3 and MCM-UA-1-0 get correlations above 0.2. The standardized deviations for vertical pumping have a quite wide range, and most of the models simulate a stronger amplitude than that of Argo data. In particularly, the amplitude of CIESM is 4.5 times that of the reference data. Nevertheless, EC-Earth3-Veg and FGOALS-g3 simulate a similar amplitude of vertical pumping with that in Argo. For the lateral induction term, the correlations vary from 0.19 to 0.47. MME outperforms all single model, with a correlation of 0.47. Ten models including CESM2-FV2, CIESM, EC-Earth3-Veg, FGOALS-g3, FIO-ESM-2-0, GISS-E2-1-G, IPSL-CM6A-LR, NESM3, NorESM2-LM and SAM0-UNICON get correlations higher than 0.30, accounting for nearly half of the analyzed models. The standardized deviations for lateral induction are also more reasonable than that of vertical pumping term, with almost models located in the range of 0.5−1.5. CESM2, FGOALS-g3, and IPSL-CM6A-LR get better scores in this term. The situation with reference to SODA is similar to that referring to Argo, the simulation skill for lateral induction term is better both in correlation coefficient and standardized deviation.
Figure
11.
Taylor diagrams for vertical pumping and lateral induction simulations for the 21 Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME) with Argo (a) and Simple Ocean Data Assimilation (SODA) (b) data.
According to the analysis mentioned above, the excellent models are different for the simulations of vertical pumping term, lateral induction term and their ratios. The reasons for the good performance of GISS-E2-1-G in simulating subduction rate mentioned above are analyzed preliminarily. GISS-E2-1-G gets high score in subduction rate, because the simulation skills are all in the forefront regardless of the MLD, vertical pumping, lateral induction or the ratio between the two compared with the reference data. In addition, FGOALS-g3 also performs well when compared with Argo, but the ratio between lateral induction and vertical pumping is larger than that of GISS-E2-1-G.
Compared with the subduction rate simulated in CMIP5, the subduction rate centers are presented better in the CMIP6, regarding the position and its two components, especially in the Southern Ocean. The mode water is generally poorly simulated in CMIP5 models. Low potential vorticity parcels and very different places in the water column are found in the Southern Ocean compared with the real ocean. In addition, the models subduct water correctly through Ekman pumping, but underestimate the larger volumes of water through the lateral induction term (Sallée et al., 2013a). However, in the CMIP6 models, the simulations are improved. Moreover, the lateral induction skills are better than that of vertical pumping with reference to Argo (Fig. 11a).
4.
Long-term trend
4.1
Time variation trend of the subduction rate
The linear trend of the global subduction rate is analyzed in detail. Due to the short time span of Argo, SODA and ORA-S3 are employed here for the long term trend analysis. Figure 12 shows the global average subduction rate anomaly of the SODA data, ORA-S3 data, 21 CMIP6 models and MME. The climatology periods are 1901−2009, 1959−2010, and 1901−2013 for the SODA data, ORA-S3 data and CMIP6 models, respectively. There is an obvious ascending trend in the global average subduction rate according to the SODA data, with an upward rate of 7.31 m/a every 100 a. From the interdecadal scale, before the 1960s, the global average subduction rate remained a negative anomaly. After that, the anomaly changed from negative to positive, with an obvious upward trend. As known well, SODA is reanalysis dataset. The long-term trend of subduction rate reflected in SODA needs additional evidence from other datasets for mutual confirmation. The trend in ORA-S3 is also analyzed here. The results show that the ascending trend is robust in both SODA and ORA-S3 data, the slopes are 13.27 m/a and 11.56 m/a per 100 a at the period of 1959−2010, respectively. Moreover, the long-term change characteristic has been confirmed by other studies. For example, Liu and Wu (2012) reached a similar conclusion. Unfortunately, none of the 21 CMIP6 models can simulate the upward trend of the subduction rate as strong as that derived from the SODA data. Among them, the long-term variation trend of the subduction rate simulated by IPSL-CM6A-LR is positive with a slope of 1.19 m/a per 100 a. The linear trends of the subduction rate simulated by the other models are negative. The time variation trend of the MME is −0.66 m/a per 100 a. There are differences in the linear trend among these models. These values vary from −4.85 m/a to 1.19 m/a every 100 a, which are very weak compared with the trend of SODA data.
Figure
12.
The global mean subduction rate time anomaly of the Simple Ocean Data Assimilation (SODA) data, ORA-S3 data, Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME).
To further clarify the spatial distribution of the long-term linear trend of the global subduction rate, we checked the long-term trend distributions of the SODA data and model simulations (Figs 13 and 14). According to the SODA data, the upward and downward trends in the 20th century coexisted, and upward trends predominated. The ascending trend that is significant at the 95% level mainly distributes in the high centers of subduction, including the northern area of the North Pacific and the North Atlantic, as well as the sea area near the ACC. The ACC area is the main region contributing to the linear positive variation trend of global subduction. The descending trend of the subduction rate is mainly distributed in the North Indian Ocean, the central North Pacific and parts of the North Atlantic, and only sporadic areas pass the significance test. In the MME results, the negative trend accounts for a large proportion, and the significant trends of the increasing subduction rate are not reflected in the northern North Pacific, northern North Atlantic or ACC region.
Figure
13.
Spatial distribution of the linear trend of the subduction rate from the Simple Ocean Data Assimilation (SODA) data (a) and multi-model ensemble mean (MME) (b). The black dots indicate that the area passed the 95% significance test.
Figure
14.
The subduction rate time trend distributions of the Coupled Model Intercomparison Project Phase 6 (CMIP6) models. The black dots indicate that the area passed the 95% significance test.
Figure 14 shows the spatial distribution of the linear trend of the subduction rate for the 21 CMIP6 climate models. Compared with the SODA data, a common problem is that the linear trends of the CMIP6 models are relatively small, and only sporadic areas in most models pass the significance test. According to the significant key sea areas in the SODA data, EC-Earth3-Veg, FGOALS-f3-L, FGOALS-g3, IPSL-CM6A-LR, MPI-ESM1-2-LR, NESM3 and NorESM2-LM can reproduce at least one of the key regions and perform better in depicting the linear trend distribution of the subduction rate compared with the other models. Among them, IPSL-CM6A-LR has a large area of significant positive trend in the southern Indian Ocean, and the positive trends in the northern North Atlantic are also reflected, which is closer to the corresponding distribution in the SODA data.
By analyzing the linear trends of the SODA data and CMIP6 models, it can be found that the largest positive trends of the global subduction rate are mainly distributed in the ACC region. Therefore, whether the variation trends in this region can be well simulated is critical for the model.
4.2
Long-term trends in the subduction rate for the ACC region
Since variation trends in subduction in the ACC region (35°−65°S, as shown in Fig. 15a) are an important piece to distinguish the model’s ability to simulate long-term trends in the global subduction rate, we further analyzed temporal changes in this region.
Figure
15.
The latitude bands analyzed in the Antarctic Circumpolar Current (ACC) region (a), the black bold lines indicate the boundaries; time variation curves of the subduction rate, vertical pumping term and lateral induction term in the ACC region for the Simple Ocean Data Assimilation (SODA) data (b).
By analyzing temporal changes in the subduction rate, the vertical pumping term and the lateral induction term of the SODA data in the ACC region (Fig. 15b), it can be seen that regardless of the long-term trend or time variation amplitude, the changes in the subduction rate are mainly dominated by the lateral induction term in this region, and the vertical pumping term shows little influence.
Next, the lateral induction term simulated by the CMIP6 models is evaluated. Figure 16 shows the time variation curves of the lateral induction term for the SODA data and CMIP6 models in ACC region. Similar to the subduction rate, most of the models are not good at simulating the ascending trend of the lateral induction term. The trends are positive for EC-Earth-Veg, IPSL-CM6A-LR and NESM3. Among them, IPSL-CM6A-LR (3.3 m/a every 100 a) is the model with the largest trend values which is an order of magnitude larger than that in EC-Earth-Veg and NESM3.
Figure
16.
Time variation of the lateral induction term (TVLI) for the Simple Ocean Data Assimilation (SODA) data and Coupled Model Intercomparison Project Phase 6 (CMIP6) models. The solid line shows the SODA data results, the dotted line shows the model simulation results, and the value in the upper left corner is the linear slope of the model result.
The lateral induction term of the subduction rate is mainly affected by the horizontal velocity at the BML and the horizontal gradient of the MLD, which represents the contribution to the subduction process due to the inclination of the mixed layer under the effect of the background flow field. To determine the simulation problem for the CMIP6 models on the lateral induction term, the horizontal velocity at the BML and the horizontal gradient of the MLD in the ACC region are analyzed in detail. As subduction mainly occurs in winter, September is chosen here.
Figure 17 shows the average zonal velocity and the average meridional velocity at the BML, and the time variations of the zonal gradient and the meridional gradient of the MLD in the ACC region from the SODA data, CMIP6 models and MME in September. Among these variables, the time changes of the meridional velocity at the BML and the meridional gradient of the MLD are relatively stable (Figs 17b, d), and the MME results are in good agreement with the SODA data. The variation amplitude of the individual models is larger than that of the SODA data, which is manifested in the meridional velocity amplitude at the BML of MCM-UA-1-0 being obviously larger. The meridional gradient of the MLD of FGOALS-g3 and NESM3 is obviously greater than that of the SODA data. The difference between each model and the SODA data is mainly reflected in the zonal velocity and zonal gradient of the mixed layer. As shown in Fig. 17a, similar to the time point of decadal changes in the global subduction rate, the decline for the zonal velocity at the BML in the ACC region appeared in the mid-1960s and accelerated in the 1990s. The zonal gradient of the MLD in this area also showed a downward trend in the corresponding time period but with a less obvious trend. Therefore, the CMIP6 models and MME show poor ability in simulating the interdecadal variation trends in the zonal velocity at the BML and zonal gradient of the MLD, which should be the main reason why the CMIP6 models cannot simulate the variation trends in the regional subduction rate of the ACC region.
Figure
17.
Abnormal time variation curves of the zonal velocity at the bottom mixed layer (BML) (a), the meridional velocity at the BML (b), the zonal gradient of the mixed layer depth (MLD) (c) and the meridional gradient of the MLD for the Simple Ocean Data Assimilation (SODA) data, Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME) in the Antarctic Circumpolar Current (ACC) region during September (d).
We comprehensively evaluated the performance of 21 CMIP6 coupled global climate models in simulating subduction rate climatology and their capacities in reproduction for the MLD, vertical pumping and lateral induction, performed trend analyses of the subduction rate and compared them with Argo and SODA data. In terms of the climatology simulation, GISS-E2-1-G can get high scores in all aspects, regardless of MLD, subduction rate, vertical pumping, lateral induction or the ratio between the lateral induction and vertical pumping. The ocean component in GISS-E2-1-G is GISS Ocean V1. Kelley et al. (2020) noted that the ocean model of GISS-E2-1-G is greatly improved in terms of the mesoscale vortex parameterization, vertical mixing, advection scheme and vertical resolution. All these improvements are closely associated with the subduction process, which may contribute to the good performance in this regard. In terms of the long-term trend of the subduction rate, none of the models can reproduce the positive long-term trend of the global subduction rate as strong as that in SODA data. IPSL-CM6A-LR simulates a positive trend with a slope of 1.19 m/a per 100 a and the slopes in other models are all negative. Through further analysis of the time variation trend of the subduction rate in the ACC region and the corresponding lateral induction term, it is still IPSL-CM6A-LR who simulates the largest trend for the long-term trend of the subduction rate in the ACC region.
In terms of the long-term variation trends of the global subduction rate, our analysis shows that the ACC region is a key area for the positive long-term trend of the global subduction rate, and in this region, subduction is mainly dominated by the lateral induction term, while the influence of the vertical pumping term can be ignored. The lateral induction represents the contribution due to the inclination of the mixed layer and the background flow field. The analysis of the time variation of the horizontal velocity at the BML and the horizontal gradient of the MLD in the ACC region during September shows that the poor simulation ability of the CMIP6 models for the zonal velocity at the BML and the zonal gradient of the MLD in this region is the reason why the models cannot reproduce the strong positive trend of the subduction rate.
In addition, a number of recent studies have pointed out that mesoscale to submesoscale processes in the ocean have important effects on mode water and subduction processes (Xu et al., 2014, 2017; Liu et al., 2019). Xu et al. (2016) studied the effects of ocean mesoscale eddies on the North Pacific Subtropical Mode water on the south side of the Kuroshio Extension based on buoy observations. Their study suggests that ocean eddies produce low potential vorticity water transport across isopycnals, causing the southward subduction of this water body into the thermocline and mode water formation. Furthermore, the contribution of these eddies to the subduction of low potential vorticity water in the mixed layer is quite considerable, accounting for more than half of its total subduction. Herraiz-Borreguero and Rintoul (2010) pointed out that advection of anomalous water from eddies and meanders of fronts bounding the Subantarctic Zone can influence the SAMW. The ACC is one of the strongest current systems on Earth, and the unstable shear of its mean flow is very conducive to the generation of mesoscale eddies. Therefore, we conclude that mesoscale and submesocale processes have important effects on the subduction of this region. However, limited by computing resources and spatial resolution, the CMIP6 models are not yet able to simulate the mesoscale to submesoscale processes, which may be one direction to improve the subduction simulation capabilities of these models in the future.
Under the global warming background, it is indicated that both observed and simulated vertical stratification in the upper ocean would be increased (Li et al., 2020; Sallée et al., 2013b). Stable stratification inhibits vertical mixing and impacts the vertical exchanges of heat, carbon, dissolved oxygen and nutrients. The MLD that is governed by the balance of stratification and buoyancy flux might be shallower in the future (Qiu et al., 2021). The volume and properties of mode waters are tightly linked to the mixed layer properties. It could be expected that the enhanced stratification and the associated decrease in ocean mixing layer will be consistent with a decline in the subduction rate eventually.
Bates N R, Pequignet A C, Johnson R J, et al. 2002. A short-term sink for atmospheric CO2 in subtropical mode water of the North Atlantic Ocean. Nature, 420(6915): 489–493. doi: 10.1038/nature01253
Chen Xingrong, Liu Shan, Cao Yi, et al. 2018. Potential effects of subduction rate in the key ocean on global warming hiatus. Acta Oceanologica Sinica, 37(3): 63–68. doi: 10.1007/s13131-017-1130-z
Chen Ju, Qu Tangdong, Sasaki Y N, et al. 2010. Anti-correlated variability in subduction rate of the western and eastern North Pacific Oceans identified by an eddy-resolving ocean GCM. Geophysical Research Letters, 37(23): L23608. doi: 10.1029/2010GL045239
Eyring V, Bony S, Meehl G A, et al. 2016. Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization. Geoscientific Model Development, 9(5): 1937–1958. doi: 10.5194/gmd-9-1937-2016
Gao Libao, Rintoul S R, Yu Weidong. 2018. Recent wind-driven change in subantarctic mode water and its impact on ocean heat storage. Nature Climate Change, 8: 58–63. doi: 10.1038/s41558-017-0022-8
Giese B S, Ray S. 2011. El Niño variability in simple ocean data assimilation (SODA), 1871–2008. Journal of Geophysical Research: Oceans, 116(C2): C02024. doi: 10.1029/2010JC006695
Gu Daifang, Philander S G H. 1997. Interdecadal climate fluctuations that depend on exchanges between the Tropics and Extratropics. Science, 275(5301): 805–807. doi: 10.1126/science.275.5301.805
Herraiz-Borreguero L, Rintoul S R. 2010. Subantarctic mode water variability influenced by mesoscale eddies south of Tasmania. Journal of Geophysical Research: Oceans, 115(C4): C04004. doi: 10.1029/2008JC005146
Hong Yu, Du Yan, Xia Xingyue, et al. 2021. Subantarctic mode water and its long-term change in CMIP6 models. Journal of Climate, 34(23): 9385–9400. doi: 10.1175/JCLI-D-21-0133.1
Huang Ruixin, Qiu Bo. 1994. Three-dimensional structure of the wind-driven circulation in the subtropical North Pacific. Journal of Physical Oceanography, 24(7): 1608–1622. doi: 10.1175/1520-0485(1994)024<1608:TDSOTW>2.0.CO;2
Kelley M, Schmidt G A, Nazarenko L S, et al. 2020. GISS-E2.1: configurations and climatology. Journal of Advances in Modeling Earth Systems, 12(8): e2019MS002025. doi: 10.1029/2019MS002025
Kubokawa A. 1999. Ventilated thermocline strongly affected by a deep mixed layer: a theory for subtropical countercurrent. Journal of Physical Oceanography, 29(6): 1314–1333. doi: 10.1175/1520-0485(1999)029<1314:VTSABA>2.0.CO;2
Ladd C, Thompson L A. 2002. Decadal variability of North Pacific Central Mode Water. Journal of Physical Oceanography, 32(10): 2870–2881. doi: 10.1175/1520-0485(2002)0322.0.CO;2
Levitus S. 1982. Climatological Atlas of the World Ocean. Princeton, NJ: NOAA
Li Guancheng, Cheng Lijing, Zhu Jiang, et al. 2020. Increasing ocean stratification over the past half-century. Nature Climate Change, 10(12): 1116–1123. doi: 10.1038/s41558-020-00918-2
Liu Zhengyu, Huang Boyin. 1998. Why is there a tritium maximum in the central equatorial Pacific thermocline?. Journal of Physical Oceanography, 28(7): 1527–1533. doi: 10.1175/1520-0485(1998)028<1527:WITATM>2.0.CO;2
Liu Lingling, Huang Ruixin. 2012. The global subduction/obduction rates: Their interannual and decadal variability. Journal of Climate, 25(4): 1096–1115. doi: 10.1175/2011JCLI4228.1
Liu Lingling, Wang Fan, Huang Ruixin. 2011. Enhancement of subduction/obduction due to hurricane-induced mixed layer deepening. Deep-Sea Research Part I: Oceanographic Research Papers, 58(6): 658–667. doi: 10.1016/j.dsr.2011.04.003
Liu Chengyan, Wu Lixin. 2012. An intensification trend of South Pacific Mode Water subduction rates over the 20th century. Journal of Geophysical Research: Oceans, 117(C7): C07009. doi: 10.1029/2011JC007755
Liu Cong, Xu Lixiao, Xie Shangping, et al. 2019. Effects of anticyclonic eddies on the multicore structure of the North Pacific subtropical mode water based on Argo observations. Journal of Geophysical Research: Oceans, 124(11): 8400–8413. doi: 10.1029/2019JC015631
Luo Yiyong, Liu Qinyu, Rothstein L M. 2009. Simulated response of North Pacific mode waters to global warming. Geophysical Research Letters, 36(23): L23609. doi: 10.1029/2009GL040906
Luo Yiyong, Liu Qinyu, Rothstein L M. 2011. Increase of South Pacific eastern subtropical mode water under global warming. Geophysical Research Letters, 38(1): L01601. doi: 10.1029/2010GL045878
Ma Jie, Lan Jian. 2017. Interannual variability of Indian Ocean subtropical mode water subduction rate. Climate Dynamics, 48(11): 4093–4107. doi: 10.1007/s00382-016-3322-1
McPhaden M J, Zhang Dongxiao. 2002. Slowdown of the meridional overturning circulation in the upper Pacific Ocean. Nature, 415(6872): 603–608. doi: 10.1038/415603a
Oka E, Qiu Bo. 2012. Progress of North Pacific mode water research in the past decade. Journal of Oceanography, 68(1): 5–20. doi: 10.1007/s10872-011-0032-5
Oka E, Suga T. 2005. Differential formation and circulation of North Pacific central mode water. Journal of Physical Oceanography, 35(11): 1997–2011. doi: 10.1175/JPO2811.1
Palter J B, Lozier M S, Barber R T. 2005. The effect of advection on the nutrient reservoir in the North Atlantic subtropical gyre. Nature, 437(7059): 687–692. doi: 10.1038/nature03969
Qiu Zishan, Wei Zexun, Nie Xunwei, et al. 2021. Southeast Indian Subantarctic Mode water in the CMIP6 coupled models. Journal of Geophysical Research: Oceans, 126(7): e2020JC016872. doi: 10.1029/2020JC016872
Qu Tangdong, Chen Ju. 2009. A North Pacific decadal variability in subduction rate. Geophysical Research Letters, 36(22): L22602. doi: 10.1029/2009GL040914
Qu Tangdong, Xie Shangping, Mitsudera H, et al. 2002. Subduction of the North Pacific mode waters in a global high-resolution GCM. Journal of Physical Oceanography, 32(3): 746–763. doi: 10.1175/1520-0485(2002)032<0746:SOTNPM>2.0.CO;2
Sallée J B, Shuckburgh E, Bruneau N, et al. 2013a. Assessment of Southern Ocean water mass circulation and characteristics in CMIP5 models: Historical bias and forcing response. Journal of Geophysical Research: Oceans, 118(4): 1830–1844. doi: 10.1002/jgrc.20135
Sallée J B, Shuckburgh E, Bruneau N, et al. 2013b. Assessment of Southern Ocean mixed-layer depths in CMIP5 models: Historical bias and forcing response. Journal of Geophysical Research: Oceans, 118(4): 1845–1862. doi: 10.1002/jgrc.20157
Suga T, Aoki Y, Saito H, et al. 2008. Ventilation of the North Pacific subtropical pycnocline and mode water formation. Progress in Oceanography, 77(4): 285–297. doi: 10.1016/j.pocean.2006.12.005
Toyama K, Iwasaki A, Suga T. 2015. Interannual variation of annual subduction rate in the North Pacific estimated from a gridded Argo product. Journal of Physical Oceanography, 45(9): 2276–2293. doi: 10.1175/JPO-D-14-0223.1
Xia Xingxue, Xu Lixiao, Xie Shangping, et al. 2021. Fast and slow responses of the Subantarctic Mode Water in the South Indian Ocean to global warming in CMIP5 extended RCP4.5 simulations. Climate Dynamics, 56(9): 3157–3171. doi: 10.1007/s00382-021-05635-w
Xu Lixiao, Li Peiliang, Xie Shangping, et al. 2016. Observing mesoscale eddy effects on mode-water subduction and transport in the North Pacific. Nature Communications, 7: 10505. doi: 10.1038/ncomms10505
Xu Lixiao, Xie Shangping, Jing Zhao, et al. 2017. Observing subsurface changes of two anticyclonic eddies passing over the Izu-Ogasawara Ridge. Geophysical Research Letters, 44(4): 1857–1865. doi: 10.1002/2016GL072163
Xu Lixiao, Xie Shangping, McClean J L, et al. 2014. Mesoscale eddy effects on the subduction of North Pacific mode waters. Journal of Geophysical Research: Oceans, 119(8): 4867–4886. doi: 10.1002/2014JC009861
Figure 1. The mixed layer depth for Argo in spring (a), summer (b), autumn (c) and winter (d); e−h are the same as a−d, but for Simple Ocean Data Assimilation results; i−l are the results for multi-model ensemble mean.
Figure 2. Taylor diagram for mixed layer depth in four seasons for the 21 Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME) with Argo (a) and Simple Ocean Data Assimilation (SODA) (b) data.
Figure 3. Zonal averaged mixed layer depth (MLD) during the spring (a), summer (b), autumn (c) and winter (d) for the Coupled Model Intercomparison Project Phase 6 (CMIP6) models, Simple Ocean Data Assimilation (SODA) data, multi-model ensemble mean (MME), as well as the reference datasets.
Figure 4. Subduction rate climatology of Argo data (a) and Simple Ocean Data Assimilation (SODA) data (b).
Figure 5. Subduction rate climatology for each Coupled Model Intercomparison Project Phase 6 (CMIP6) model.
Figure 6. Taylor diagram of subduction rate climatology simulation of the 21 Coupled Model Intercomparison Project Phase 6 (CMIP6) models and the multi-model ensemble mean (MME) with Argo (a) and Simple Ocean Data Assimilation (SODA) (b) data.
Figure 7. Subduction rate climatology of the Coupled Model Intercomparison Project Phase 6 (CMIP6) multi-model ensemble mean (MME) (a), the bias relative to Argo (b) and the bias relative to Simple Ocean Data Assimilation (SODA) (c).
Figure 8. Vertical pumping (a) and lateral induction (b) climatology of Argo; c and d are the corresponding results of Simple Ocean Data Assimilation (SODA) data. The value in the top right-hand corner represents its contribution to subduction rate.
Figure 9. Vertical pumping climatology for each Coupled Model Intercomparison Project Phase 6 (CMIP6) model. The value in the top right-hand corner represents its corresponding contribution to subduction rate.
Figure 10. Lateral induction climatology for each Coupled Model Intercomparison Project Phase 6 (CMIP6) model. The value in the top right-hand corner represents its corresponding contribution to subduction rate.
Figure 11. Taylor diagrams for vertical pumping and lateral induction simulations for the 21 Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME) with Argo (a) and Simple Ocean Data Assimilation (SODA) (b) data.
Figure 12. The global mean subduction rate time anomaly of the Simple Ocean Data Assimilation (SODA) data, ORA-S3 data, Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME).
Figure 13. Spatial distribution of the linear trend of the subduction rate from the Simple Ocean Data Assimilation (SODA) data (a) and multi-model ensemble mean (MME) (b). The black dots indicate that the area passed the 95% significance test.
Figure 14. The subduction rate time trend distributions of the Coupled Model Intercomparison Project Phase 6 (CMIP6) models. The black dots indicate that the area passed the 95% significance test.
Figure 15. The latitude bands analyzed in the Antarctic Circumpolar Current (ACC) region (a), the black bold lines indicate the boundaries; time variation curves of the subduction rate, vertical pumping term and lateral induction term in the ACC region for the Simple Ocean Data Assimilation (SODA) data (b).
Figure 16. Time variation of the lateral induction term (TVLI) for the Simple Ocean Data Assimilation (SODA) data and Coupled Model Intercomparison Project Phase 6 (CMIP6) models. The solid line shows the SODA data results, the dotted line shows the model simulation results, and the value in the upper left corner is the linear slope of the model result.
Figure 17. Abnormal time variation curves of the zonal velocity at the bottom mixed layer (BML) (a), the meridional velocity at the BML (b), the zonal gradient of the mixed layer depth (MLD) (c) and the meridional gradient of the MLD for the Simple Ocean Data Assimilation (SODA) data, Coupled Model Intercomparison Project Phase 6 (CMIP6) models and multi-model ensemble mean (MME) in the Antarctic Circumpolar Current (ACC) region during September (d).
DownLoad:
DownLoad:
DownLoad:
