1 INTRODUCTION

An important part of the ventilated thermocline theory (e.g., Iselin, 1939; Stommel, 1979; Luyten et al., 1983; Pedlosky and Young, 1983; Huang, 1990) is the subduction process of upper ocean water. The subduction process refers to the ocean dynamic process in which the water in the mixed layer is detrained into the permanent thermocline along the isopycnal surface under the effect of buoyancy forcing, for example, surface cooling, and flows into the subsurface and deeper layers to participate in thermohaline circulation. Because the abrupt deepening of the mixed layer mainly occurs in late winter, subduction is taken as an intermittent process (Stommel, 1979; Williams et al., 1995), in which subducted water is selected by a "demon", which is known as "Stommel's demon". Cushman-Roisin (1987) further divided the water exchange between the mixed layer and thermocline into three periods: late winter to early spring, early spring to early autumn, and early autumn to late winter. During the period from early spring to early autumn, the water detrained from the mixed layer usually only enters the seasonal thermocline. When the mixed layer becomes deeper, this water is most likely re-entrained back into the mixed layer; therefore, this is not effective subduction. During the period from early autumn to late winter, along with the deepening of the mixed layer, the subsurface water continuously entrains into the mixed layer through the bottom of the mixed layer. This process is called entrainment. According to Stommel's hypothesis, the short period between late winter and early spring is a period of effective subduction when the water in the mixed layer enters the permanent thermocline where is too deep for water to re-entrain into the mixed layer. Since only the effectively subducted water mass will carry the surface seawater characteristics into the subsurface layer and then into the meridional overturning circulation (MOC), it is particularly important to estimate the effective subduction rate precisely.

In general, the mode water originating from the subduction process is considered the source of ocean thermohaline circulation and has a very important contribution to its formation and variation. Therefore, many previous studies have focused on the subduction process of mode waters. For example, in the subtropical gyre, subducted mode waters transport heat and salt from the subtropics to the tropics (McCreary and Lu, 1994). In the subtropical North Atlantic, the large surface density due to the high salinity leads to strong subduction and forms the Subtropical Mode Water (STMW), which thereby participates in the MOC of the North Atlantic (e.g., Häkkinen, 1999; Schott et al., 2003; Zhang et al., 2003) and regulates the global climate. In the subtropical Pacific, the subduction rate of STMWs and Tropical Waters (TWs) and their transport to the equatorial region affect the interannual variability in the El Niño-Southern Oscillation (ENSO) and even the global climate (e.g., Kleeman et al., 1999; Qu et al., 2008; Liu and Wu, 2012; Nie et al., 2016; Tsubouchi et al., 2016; Qu and Gao, 2017). Therefore, the evaluation of the effective subduction rate in different regions is an important approach to study the formation of different mode waters and their variability, and it is also an important way for us to obtain a deeper understanding of the global climate system.

Since the mixed layer in most regions has a significant seasonal cycle, the estimation of the effective subduction rate in most regions is usually calculated in the form of the annual subduction rate (Qiu and Huang, 1995). There are two main methods for calculating the annual subduction rate: one is in the Eulerian frame (henceforth referred to as the "Eulerian method"), and the other is in the Lagrangian frame. The Eulerian method first determines the period of effective subduction within a year and then integrates the local subducted water in the effective period to obtain the annual subduction rate (e.g., De Szoeke, 1980; Cushman-Roisin, 1987; Da Costa et al., 2005; Liu et al., 2016; please see Section 2 for details). This method is accurate in theory; however, it has a great limitation in practical applications: the data used must have a high time resolution. If monthly data are adopted, the accuracy of the effective subduction period is only one month. Most subduction processes occur within 1–2 months during the deep winter (Stommel, 1979; Williams et al., 1995); therefore, the effective period obtained by monthly data is inaccurate. At present, the best observation data that can be used for ocean subduction analysis are gridded Argo data. Since the temporal resolution of Argo data is monthly, Argo data are obviously unsuitable for use in the Eulerian method.

Therefore, the practical method to calculate the annual subduction rate using observation data is the Lagrangian particle tracking method (henceforth referred to as the "Lagrangian method"), as in Qiu and Huang (1995). The Lagrangian method equation for calculating the annual subduction rate S_ann is as follows:

The first two terms on the right are vertical pumping (VP), which represents the effect of the vertical velocity w_mb at the bottom of the mixed layer. The volume of water detrained from the bottom of the mixed layer to the thermocline can then be estimated with w_mb; w_Ek is the Ekman pumping term; and is the correction term by the β spiral, which is the difference between the velocity at the bottom of the Ekman layer and that at the bottom of the mixed layer. β is the gradient in the planetary vorticity, and the overbar of the term refers to the average along the Lagrangian trajectory of particles over one year. The third term is lateral induction (LI), in which h_{m, 0} and h_{m, 1} are the mixed layer depths (MLDs) of the first winter and the second winter, respectively, and T is 1 year. The physical significance of the LI term is that the particles at the bottom of the mixed layer will flow downstream by advection. If the MLD downstream is relatively shallow, the particles will reach a deeper position than the bottom of the mixed layer, which means that there will be more water detrained at the original position. The difference between these two MLDs is the subduction component caused by advection. Due to its clear theory and convenient calculation, the Lagrangian method is now widely used by oceanographers (e.g., Qiu and Huang, 1995; Karstensen and Quadfasel, 2002; Qu et al., 2008, 2017; Liu and Huang, 2012; Liu and Wu, 2012; Mill and De Moraes Paiva, 2013).

However, the Lagrangian method has its own problems. In addition to ignoring the mixing effect (Gao et al., 2011), there are still two noteworthy problems:

1. It assumes that the volume of effective subduction at a given source location is determined by the depth of the particles released at the bottom of the MLD in winter at the source location, which could finally reach a deeper position after one year of tracking, but the actual subduction volume at the source location is not considered.

2. It assumes that Stommel's demon works everywhere (Stommel, 1979). That is, subduction only occurs during a very short period in which the mixed layer rapidly becomes shallow in the deep winter. Therefore, we consider that the water column is almost simultaneously subducted at the source and maintains a straight translational motion at all times (Fig. 1).

Corresponding to the two problems, this method may have two deviations:

1. In fact, the MLD has high spatial variability. If the variation in the MLD at the source location is relatively small, the corresponding possible subduction should also be relatively small. Furthermore, if the winter MLD at the terminal location is very shallow, or the downward Ekman pumping along the particle trajectory is very strong, or the effects of both situations exist at the same time, the subduction volume calculated by the Lagrangian method is likely to be greater than that actually subducted at the source, which will lead to the overestimation of the subduction rate (see Section 3.1 for details).

2. According to a previous study (Williams et al., 1995), Stommel's demon works in the subtropical gyre. However, according to the definition of subduction, effective subduction may not occur only in the deep winter. Imagine some regions where the wind is sufficiently strong, leading to a large downward Ekman pumping; therefore, effective subduction is also likely to occur in periods other than deep winter (see Section 3.2 for details), which is not considered in the Lagrangian method.

To clarify the above issues, we conducted a detailed theoretical analysis and quantitative estimation of these two possible deviations in the Southern Hemisphere. The rest of the paper is organized as follows. A brief description of the data and method of analysis is presented in Section 2. The deviations and their comprehensive deviation are theoretically described and quantitatively estimated in detail in Section 3. Results are summarized in Section 4.

2 MATERIAL AND METHOD 2.1 Argo data and NCEP/NCAR reanalysis wind product

Argo data are one of the most important types of ocean observation data at present. Thousands of Argo floats have been deployed since the early 2 000 s to obtain the temperature and salinity (T/S) of the ocean from a typical upper level of 5 m to 2 000 m. By using all qualified Argo profiles, a near real-time, monthly T/S product for the global ocean on a 1°×1° grid has been created by the Asian Pacific Data Research Center (APDRC) of the International Pacific Research Center, University of Hawaii. These data span from January 2005 to the present and have a total of 26 standard vertical levels consistent with those of the World Ocean Atlas (e.g., Levitus, 1983). The climatological monthly Temperature/Salinity data are obtained, from which the geostrophic current of each standard level relative to the non-motion reference level of 2 000 db is calculated. The National Center for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis wind product (Kalnay et al., 1996) is also used for the present study. By using the climatological monthly NCEP wind stress, the Ekman pumping is calculated in this study. The subduction rate is calculated by using these two datasets.

2.2 MLD

Based on the climatological Argo T/S data, the climatology of monthly MLD is calculated. The criteria used to define the MLD are diverse in previous studies (e.g., Lindstrom et al., 1987; Lukas and Lindstrom, 1991; Sprintall and Tomczak, 1992; Kara et al., 2000; De Boyer Montégut et al., 2004). According to De Boyer Montégut et al. (2004), the MLD is defined as where the potential density has increased from that of the surface by a density threshold equivalent to a temperature change of 0.5 ℃ at constant salinity. The density threshold can be expressed as

in which T_{0 m}, S_{0 m}, and P₀ are the local surface temperature, salinity, and pressure, respectively.

2.3 The jet propulsion laboratory (JPL) "off-line" simulated tracer method

The off-line simulated tracer method of Fukumori et al. (2004) is used to calculate the annual subduction rate in the Southern Hemisphere. In this method, precomputed velocity and mixing tensors are used to integrate the model's advection-diffusion equation. In the context of a general circulation model, the temporal evolution of a passive tracer is determined by the same advection-diffusion equation as that for temperature and salinity, except for differences in possible tracer sources and sinks and in external forcing and boundary conditions:

where c is the tracer concentration, u is the threedimensional velocity vector, and κ is the mixing tensor. If a particular patch of water is uniformly initialized with a passive tracer, the subsequent movement of this tracer will indicate the pathway and destination of the initial patch of water. In particular, the magnitude of the tracer at a given location relative to the value of the initial patch describes the concentration of this initial water mass at the given location. While passive tracers can trace where tracertagged water goes, the adjoint passive tracer that is integrated backward in time can identify where the adjoint tracer-tagged water comes from. Mathematically, a quantity's adjoint is the sensitivity of that quantity to other variables. In particular, the adjoint passive tracer is defined as the sensitivity of a passive tracer at a given location (target) and at a given time (terminal instant) to tracers at other locations at earlier instances. Such sensitivity will be finite if and only if water from such a location makes its way to the target at the terminal instant. Therefore, the value of this sensitivity relative to that of the target describes the fraction of water that circulates from that location to the target location. Please refer to Fukumori et al. (2004) for details.

The annual subduction rate that includes the effect of mixing in addition to advection can be analyzed using this simulated adjoint tracer method. An adjoint tracer is initialized in September (austral winter) with a unit value (in arbitrary tracer units per volume, ATU/m³) for the entire water volume under the surface mixed layer of the Southern Hemisphere (80°S–0° and 0°–360°E). Using the velocity and mixing tensors of the Estimating the Circulation and Climate of the Ocean (ECCO) simulation, the adjoint tracer is integrated backward in time for 1 year until September of the previous year. The annual subduction rate can then be calculated as

in which C(x, y, z) is the fraction of tracer-tagged water at position (x, y, z) after the 1-year backward integration, and D_m is the MLD in winter (September) of the previous year. According to this definition, if tracer-tagged water below the winter mixed layer is found re-entering the mixed layer at location A and is traced backward within the mixed layer to location B until the previous winter, it will be regarded as "subducted" at location B. Please refer to Gao et al. (2011) for details.

2.4 The Eulerian method

In this paper, the deviation in the classical Lagrangian method is corrected by using the Eulerian method. In this method, S_ann is defined by the following equation:

in which T=1 year and S_Euler is the annual subduction rate calculated by Eulerian method, is the instantaneous local detrainment rate (De Szoeke, 1980; Cushman-Roisin, 1987), where u_mb is the horizontal velocity at the bottom of the mixed layer and h_m is the instantaneous MLD. t_s and t_e are the starting and ending times of the effective detrainment, respectively, which can be determined by tracking the trajectory of particles released at every time interval from the bottom of the mixed layer and checking whether they are still under the mixed layer in the next winter. The period between t_s and t_e is the effective detrainment period.

The Eulerian method is theoretically accurate, but the problem is that the determination of its effective detrainment period depends heavily on the time resolution of the data used. Most of the effective subduction occurs during the short period of late winter to early spring, which only lasts for 1–2 months. Therefore, if only monthly data such as Argo data are used, there will be a large deviation in the determination of the effective detrainment period. Therefore, the Eulerian method can only be used to calculate the subduction rate by using model data with high time resolution.

3 RESULT

As analyzed above, two assumptions are used in the classical calculation method of the annual subduction rate:

1. The calculation of the effective subduction is determined by two factors: the new depth of particles after one-year trajectory and the MLD at the location. The amount of the subducted water at the source is not considered.

2. Subduction only occurs in the process of the mixed layer becoming shallower in the deep winter, and possible subduction due to Ekman pumping at other times is not considered.

In the following, we will perform theoretical analysis and quantitative estimation for the possible deviations caused by the above two assumptions.

3.1 Deviation caused by assumption (1)

According to Eq.1 for calculating the annual subduction rate by the classical method, a schematic diagram of the corresponding subduction process is given in Fig. 2. It is seen from the figure that the physical process of Eq.1 can be described as follows: the water particle is released at 'O' at the bottom of the mixed layer in the first winter, after moving downstream for one year the particle reaches the new position 'E'. The annual subduction rate is equal to the difference between the depth of 'E' and 'E_h', the 'E_h' is the position of MLD (h_{m, 1}) in the second winter. However, we can see an obvious problem from the figure: the height of the water column representing the subduction rate is actually larger than the real subduction at the source area: S_source. We can see from the figure that the length of the water column (O_h~O) between the initial position 'O' of the subduction and 'O_h' at the bottom of the summer MLD is the largest possible amount of subduction water at the source in winter of the first year.

When the mixed layer suddenly becomes shallower in the source area, there will be a new subducting water column, which will move horizontally due to advection until the winter of the next year and will be affected by Ekman pumping along its path, causing rise or fall in the vertical direction. As shown in Fig. 2, when the dominant direction of Ekman pumping along the path is downward, the subducting water column will sink in the process of downstream movement. However, if the downstream Ekman pumping is large or the downstream winter mixed layer is shallow, the top of the water column (E₀) may be compressed under the winter MLD when the water column reaches the position (E~E₀) deep in the second winter. Obviously, the water column (E_h~E₀) from the top of the water column to the bottom of the mixed layer is not subducted from the source area, so there will be deviation.

Therefore, in which areas will this deviation appear and is this deviation significant? A quantitative analysis of this deviation is performed in the Southern Hemisphere. It can be seen from the analysis that the amount of subducted water masses at the source is determined by three components: first, the mixed layer rapidly becomes shallower in winter, which takes place within 1–2 months after entering the deep winter. Second, the local vertical pumping at the source area also forces the water out of the mixed layer and forms subduction. Third, in these two months, due to the tilt of the mixed layer bottom, the water subducts into the permanent thermocline. Overall, the amount of subduction water masses at the source S_source can be defined as the local subduction within two months of deep winter, which can be calculated by the Eulerian method:

where t_s and t_e are the times when source subduction begins and ends, respectively. As shown in Fig. 3d, the averaged MLD of the South Pacific Ocean and the South Atlantic Ocean represented by the black dotted line begins to shallow after reaching the deepest depth in approximately September, and the period of shallowing usually takes 1 or 2 months. Although the mixed layer in many areas does not reach its shallowest level after the end of this period, the mixed layer generally shallows much more slowly in the next one or two months (Fig. 3d). We believe that after these two months, the subducted water has left the local area due to the current field. Therefore, it is reasonable to set t_s and t_e as the deepest MLD month and two months after that month, respectively; September and November are selected in most areas, but an exception occurs in the sea area south of New Zealand. Due to the influence of unique buoyancy flux characteristics, the spatial variability in MLD is obvious south of New Zealand. The red solid line with triangular symbol denoted 'A' in Fig. 3d identifies a characteristic point (53°S, 168.5°E) in the sea area south of New Zealand, which indicates that the appropriate period of rapid mixed layer shallowing should range from October to December.

Based on the above discussion, the subduction rate S_ann calculated by the classical method cannot be greater than S_source. If it is greater than the maximum subduction of the source, the excess part is the deviation in the subduction rate, that is,

where S_e is the excess part which is overestimated in the calculation of subduction rate.

In the specific analysis, first, according to Eq.1 of the classical method, the annual subduction rate and its components in the Southern Hemisphere are calculated by using the monthly averaged Argo float and NCEP wind field data. Then, according to Eqs.6 and 7, the corresponding deviation is calculated.

Figure 4 shows the spatial distribution of the climatological subduction rate and its vertical pumping and lateral induction components in the Southern Hemisphere. It can be seen from the figure that the subduction area of the South Pacific is mainly concentrated in the subpolar front (60°S–40°S) in the middle latitudes (here the range of middle latitudes is defined as 60°S–30°S) and the subtropical circulation area (25°S–10°S) in the low latitudes (the low latitude range here is 30°S–10°S). Notably, the magnitudes of subduction at different latitudes are different, and subduction at middle latitudes is much larger than that at low latitudes, which is mainly due to the large difference in MLD between low and middle latitudes (Fig. 3). From the component diagrams in Fig. 4b & c, it can be seen that the main component at middle latitudes is the lateral induction term, and vertical pumping plays a relatively small role, while the subduction at low latitudes is different, the vertical pumping term and the lateral induction term have the same magnitude, and subduction is the result of the joint influence of the two components.

Figure 5 shows the spatial distribution of S_e and the relative deviation ratio (S_e/S_ann) in the Southern Hemisphere calculated based on the above considerations. The results indicate that such a situation does exist in the ocean where the maximum subduction of the source in winter is smaller than that calculated by the classical method in some regions. It can also be seen from the contours in Fig. 5a that the significant area is mainly concentrated in the subduction area in the low latitudes, both in the Pacific and Atlantic, from approximately 22°S–10°S, and the absolute deviation in some areas is more than 50 m. The area with high deviation values in the Pacific Ocean is in the same latitude range as that in the Atlantic Ocean; the difference is that the deviation is near the east coast of the Pacific Ocean, while it is near the west coast of the Atlantic Ocean. The high deviation area mainly corresponds to the formation region of STMW characterized by high salinity in the Pacific. In the middle latitudes of the oceans, the deviation is almost zero and can be ignored. This shows that in the middle latitude area, part of the water mass subducted in winter at the source will be entrained back to the downstream mixed layer in the next winter.

To evaluate the significance of the deviation S_e, the estimated value of the relative deviation ratio S_e/S_ann in the Southern Hemisphere was calculated (Fig. 5b). It can be seen from the figure that the deviation is relatively obvious in the subduction area of the low latitudes. The relative importance in the Pacific and the Atlantic mostly exceeds 30%, and the maximum value can reach 80% (Fig. 5b). Therefore, for low latitudes, the relative deviation in the subduction rate calculated by the classical method cannot be ignored.

Additionally, the seasonal variation in the mixed layer of the Indian Ocean is totally different from that of the Pacific Ocean and the Atlantic Ocean, because it is very obvious that the Indian Ocean is affected by the monsoon, and the MLD may have more than one peak in one year (Liu et al., 2018), as shown by Fig. 3d, two peaks appear in the June (59 m) and December (48.5 m) at the represent point D marked in Fig. 3a in spite of the relative shallow MLD, in the four months after June, the mixed layer became 19 m shallower, and in the two months after December, the mixed layer became shallower by 11.5 m. The two shallowing processes of the mixed layer that occurred within one year will both contribute to the annual subduction rate. Therefore, it is not suitable to use the classical method to calculate the annual subduction rate in the middle and low latitudes of the Indian Ocean. Similarly, the above analysis method is not applicable to the middle and low latitudes of the Indian Ocean. However, a preliminary and conceptual consideration can be made: if the mixed layer has two minimum depths in one year, there will be two effective subduction periods, and the total subduction rate should be the sum of these two processes (Liu et al., 2018, 2020). In each process, the estimation method of the subduction rate is the same as the conventional method, so these two processes should both consider the deviation, and the deviation is more likely to exist in the Indian Ocean than in the Pacific and Atlantic Oceans. However, in this paper, the classical method is used to calculate the Indian Ocean subduction rate, and we do not discuss the deviation results and only give the patterns for reference.

Why does the classical method have no deviation at middle latitudes but obvious deviation at low latitudes? Through the climatological MLD during winter and summer in the Pacific Ocean (Fig. 3c), the reason can be seen. In the middle latitude area, the variation in the MLD is very intense in winter and summer. For example, in the deep winter, the maximum of the MLD at the high subduction point of the middle latitudes (101.5°W–53°S) exceeds 500 m, while in summer, the MLD is only approximately 70 m at the same location. The variation in a year is greater than 400 m, which ensures that the maximum subduction at the source area is very large. At the same time, the vertical pumping in the nearby area is relatively small (Fig. 3b), which ensures that it does not have a great impact on the vertical movement of the subducted water. When the water column reaches the downstream position not far away in the second winter, the deeper mixed layer will inevitably cut off its upper part; that is, there is no deviation Se. However, in the low latitudes, the situation is quite different. First, the difference in MLD between winter and summer is not large. For example, in the deep winter, the deepest mixed layer is 153 m at the maximum position of subduction (19°S–102.5°W) in the low latitudes, while in the summer, it is only 51 m, and the variation range is approximately 100 m (Fig. 3d). In addition, the vertical pumping in this area is much stronger than that in middle-latitude areas (Fig. 4b). As a result, the subducted water column at the source is not sufficient, and the strong vertical pumping along the path will also continuously push the subducted water column downward, and it is likely to completely compress the subducted water column below the winter mixed layer at the end of the second year; that is, deviation S_e will be generated.

The question may arise that even if the amount of water at the end exceeds that at the source, this part of the water column is under the mixed layer with certainty, so should it not belong to the effective subduction? Actually, it should be pointed out that due to the existence of the β-spiral, the horizontal current field differs strongly at different depths. Therefore, the water column (O_h~O) that subducted at the time selected by Stommel's demon should be scattered and twisted during the period of moving downstream, instead of always maintaining the original shape of the water column. Additionally, for the water column (E_h~E) at the end point calculated by the Lagrangian method, the water mass may come from the subduction at other positions except for the lowest point, but when we calculate the subduction rate, the effective subducted water mass at the initial position plays a more important role, as that area determines the temperature and salinity characteristics of the subsurface water masses. Therefore, it is not important whether the column obtained by the Lagrangian method comes from the initial position of subduction. The important feature that we care about is to see the difference between the initial subducted water column and the result of the Lagrangian method, that is, the deviation. In other words, the first deviation discussed in our paper is actually a necessary condition; that is, as long as the height of the water column (E_h~E) is greater than that of the equivalent water column (O_h~O) subducted at the source, this deviation must exist in this method.

Furthermore, the Eulerian method is directly used to estimate the water column (O_h~O), that is, the water mass subducted in the selected "Stommel's demon period" is not a cylinder but an equivalent water column height. At the same time, the span of Stommel's demon period is usually only 1–2 months. Due to the limitation of the Eulerian method in resolution, accurate effective subduction is difficult to obtain in monthly results. For example, the effective subduction period in most areas lasts approximately 1.5 months, but 1.5 months cannot be taken as the effective subduction period in the calculation. Therefore, it is impossible to use the exact effective subduction period when using the Eulerian method based on monthly averaged data such as Argo data. Therefore, when calculating the water column (O_h~O) in this study, the effective period is assumed to be two months after the deepest MLD, which is actually a maximum estimation. Therefore, if the subducted column O_h~O obtained by using 2 months as the source Stommel's demon period is less than the result obtained by the Lagrangian method (i.e., the height of the water column E_h~E), then at least the length of E_h~E is overestimated in the subduction calculation.

It should be noted that the MLD is a key factor influencing S_ann. The criteria used to define the MLD have been diverse in previous studies (e.g., Lukas and Lindstrom, 1991; Sprintall and Tomczak, 1992; Monterey and Levitus, 1997; Kara et al., 2000; De Boyer Montégut et al., 2004). In some areas where strong subduction exists, the difference of MLD defined by different criteria could be several hundred meters (Piron et al., 2016). Since the lateral induction (LI) is determined by the difference in MLD of the first and second winter, the results of LI calculated by different MLD criteria might be very different. Therefore, it is necessary to further evaluate the sensitivity of the MLD criteria in the calculation of S_ann. Due to the space constraint, it will be report in another paper.

In summary, because the classical method does not consider the source subduction, the subduction rate in some of the subduction areas in the low latitudes is overestimated, and there are two reasons for this deviation: (1) there is a large vertical pumping in the subduction area; and (2) there is a small mixed layer change at the source.

3.2 Deviation caused by assumption (2)

There is another hypothesis in the classical method: subduction only occurs in the process of the mixed layer becoming shallower in the deep winter and is not considered in other periods; this is also the main assumption of Stommel's demon. In this section, we try to modify this hypothesis based on some considerations to make it more in accordance with the actual physical process in the ocean. First, it should be noted that the vertical pumping considered in the classical method refers to the integral value (the first term of Eq.1 along the trajectory of the particle, which is used to modify the final subduction depth of the particle. However, what we should consider here is the local vertical pumping effect at the source, and the two are different. Because the main component of vertical pumping is Ekman pumping, vertical pumping is a continuous forcing, which is different from MLD variation with its sudden seasonal change. Because of the relatively large fluctuation of the mixed layer in general, the water mass subducted in any season except winter is generally considered to be re-involved in the mixed layer in the winter of the second year, i.e., invalid subduction. However, if the fluctuation in the mixed layer of a specific area is very small and the vertical pumping is very large, there is likely to be effective subduction without re-involvement.

In the following, we will give a detailed description of the mechanism by which vertical pumping contributes a part of the effective subduction rate. As shown in Fig. 6, at a source location 'O' in the Southern Hemisphere, during two months (September, October) entering the deep winter, there will be a rapid shallowing of the mixed layer (Fig. 5b). Therefore, a water mass (gray water column with dotted frame) will subduct below the mixed layer and then move downstream under the influence of advection and vertical pumping and finally reach the point 'E' in the next winter, contributing to effective subduction (solid frame gray water column). This process is the winter subduction process discussed in the previous section. However, in the following time period (November to August of the next year), there is still an opportunity for point 'O' to have effective subduction, which is not considered in the classical method. As shown in the figure, if there is strong vertical pumping locally, point O₁₁ at the bottom of the mixed layer in November may subduct below the mixed layer and move downstream along the horizontal direction. At this time, if there is still a large vertical pumping along the path, the particle will enter a deeper water layer until it reaches a new position E₁₁ in the next winter. If the MLD of this location is shallow in winter, E₁₁ will still be located below the mixed layer, which is an effective subduction. Two conditions are needed for effective subduction: (1) the local pumping and the vertical pumping along the path are relatively strong, and (2) the winter mixed layer at the end point is shallow. By analogy, if the trajectory of water particles at the bottom of the mixed layer is tracked from November to August of the next year, it can be checked whether there is effective subduction in each month. Taking the case in Fig. 6 as an example, at source point 'O', due to the large vertical pumping, water will be detrained out of the mixed layer from November to August of the next year. However, in the next winter, only the water particle from November to February is still under the winter mixed layer, so the subduction water column (thin red strip) in this period is an example of effective subduction, while the subduction water column (thin blue strip) from March to August is an example of invalid subduction.

Obviously, the analysis of the above subduction process is in accordance with the calculation method of the subduction rate in the Eulerian form, so we can use Eq.5 to estimate the quantity of subducted water in this period. Due to the complexity in each region, the months that can form such a subduction process are not necessarily continuous, as shown in Fig. 6, and there may be discontinuities. Therefore, it is necessary to judge the existence of effective subduction on a monthly basis. Our criterion is that if the water particles in a certain month and the following month are all under the winter mixed layer, then all subducted water masses in that month can be considered effective subduction. In each month, there are two local processes that determine the amount of subduction at the source point O: the local vertical pumping process and changes in the MLD. Therefore, between November and August of the next year, the deviation S_vp in the subduction rate caused by vertical pumping at the O point is equal to the sum of subducted water consisting of each effective detrained month. The discretized form of Eq.5 is:

The subscript "mon" represents the month of effective subduction except the two months that have been considered in Section 3.1. It should also be noted that since the monthly averaged data are used, the effective subduction period obtained by the Eulerian method can only be accurate to one month.

According to Eq.8, S_vp in the Southern Hemisphere is also estimated using Argo data and NCEP wind field data, and the results are shown in Fig. 7.

S_vp is mainly distributed in the low latitudes between 20°S and 10°S, which is very similar to the area of S_e in the previous section. The deviation in the east-central Pacific is the most obvious, there is a high value in the latitude band 110°W–130°W, the maximum deviation can exceed 40 m, and the maximum relative importance ratio can reach 70%. However, noticing that the locations of contours represent subduction, the main deviation is not in the typical subduction area, which indicates that the domain of the actual subduction rate is wider if we only consider this deviation. The situation in the Atlantic Ocean is similar, with the large subduction rate area located southwest of the high deviation area. The S_vp of the Indian Ocean is smaller but more widespread, and the subduction area calculated by the old method is also widespread in the Indian Ocean (Liu et al., 2018), so it is very likely to have a larger subduction rate if applying the double peak maximum MLD method. Therefore, the deviation in the case of the Indian Ocean is worthy of further investigation. Moreover, S_vp is also distributed in some areas in the middle latitudes of the Southern Ocean, especially around the sea south of New Zealand (Fig. 7a). However, it is noteworthy in Fig. 7b that the area with a high deviation ratio is not as wide as the deviation pattern, which is largely due to the high subduction rate in the high deviation area, so the influence of this deviation is not prominent. The main reason for the existence of high S_vp in the low latitudes of the Southern Ocean is that the vertical pumping in the low latitudes is strong (Fig. 4b), while the mixed layer is shallow and the fluctuation in MLD is small (Fig. 3d). Under such conditions, vertical pumping will lead to a continuous subduction process in periods other than winter. The magnitude of subduction deviation is particularly large in some low-latitude areas, where it cannot be ignored.

3.3 Comprehensive deviation and correction method

Two possible deviations in the calculation of the subduction rate by the classical method are analyzed and evaluated. As mentioned above, the conditions for the two deviations are very similar, so they are both located in the low latitudes and have similar influence areas. Therefore, if the two deviations are considered together, some of them compensate for each other, but the total deviation is still significant. We can see from Fig. 8 that the absolute deviation is mainly positive in the main subduction area, which indicates that the overestimation in the classical method of calculating the subduction rate is dominant, and the relative deviation ratio can reach 50% (Fig. 8). In the Pacific Ocean, the deviation is mainly distributed in the eastern Pacific Ocean 160°W–90°W, 20°S–10°S, which corresponds to the subduction area of Eastern Subtropical Mode Water (ESTMW) in the eastern South Pacific. Previous studies have shown that the subduction of ESTMW and its transport to the equator will have an important modulation effect on the climate of the tropical region (Liu and Wu, 2012; Nie et al., 2016; Qu and Gao, 2017), so the deviation in the calculation of the annual subduction rate in this region using the classical method will affect our understanding of this process (Liu and Wu, 2012). In the low latitudes of the Atlantic Ocean, the overestimation of the subduction rate is also obvious, and high values of (S_e–S_vp)/S_ann appear in the area close to the equator. This result reminds us that when using the classical method to calculate the subduction rate, the subduction area in the low latitudes may be inaccurate.

To further test the above viewpoint, we use the ECCO model data from the Jet Propulsion Laboratory (JPL) to calculate the subduction rates in the South Pacific Ocean and the South Atlantic Ocean through the classical method and the simulated passive tracer method and compare the two results (Fig. 9). Because the classical method does not consider the role of mixing within the theoretical framework, here, we only show the no-mixing result of the latter method for comparison; more details can be found in Gao et al. (2011). The results clearly show that the two patterns (Fig. 9a & b) are very similar, indicating that the calculation of the subduction rate using the classical method is simple and accurate. However, in the low latitudes of the South Pacific and South Atlantic, the subduction rate obtained by the classical method is obviously larger than that obtained by the simulated passive tracer method, which is consistent with the above deviation analysis (Fig. 9c). Since the subduction rate obtained by the simulated passive tracer method can be regarded as an accurate result, the analysis results of the above two deviations in this paper can explain the difference very well.

In summary, according to the above analysis and estimation, we can use the Eulerian method to make a preliminary correction for the classical method. The steps are summarized as follows:

1. Use Eq.1 of the classical method to calculate the annual subduction rate S_ann.

2. Calculate the local maximum detrained water volume in two months entering the deep winter by using the Eulerian method Eq.6.

3. Use Eq.8 of the Eulerian method to calculate the effective subduction S_vp in the following months.

4. Use the two results obtained in steps 2 and 3 to modify the annual subduction rate of the classical method S_ann:

Please note that the Lagrangian method is not used in Stommel's demon period in the first case, so in areas where both deviations exist, the revised Equation is equivalent to the Eulerian calculation method of the subduction rate: S_Euler.

4 DISCUSSION

As a process of ventilation, subduction relates the properties of the atmosphere with the subsurface water properties of the ocean, which has a considerable impact on the formation and transport of mode water in many regions. Among the world oceans, the Antarctic Circumpolar Current (ACC) plays an important role as a conduit for the interbasin exchange of water masses (Weijer et al., 2012). Through the concept of the Southern Hemisphere super gyre (SHSG), the water properties of the whole Southern Ocean also link together. Then, by the MOC of each ocean, especially the MOC of the Atlantic Ocean, the water masses are globally redistributed. The major water masses involved in these processes, such as Subantarctic Mode Water (SAMW), Antarctic Intermediate Water (AAIW), and Antarctic Bottom Water (AABW), are of great importance to climate change. For example, Gao et al. (2018) pointed out that the heat content related to the global warming hiatus is stored in SAMW and AAIW. Most previous studies focused on the inherent properties of water masses, such as the thickness, depth, and heat content. By calculating the subduction rate, there will be a new way to study the relationship between water masses and climate change from another perspective. Downes et al. (2010) pointed out that the volume of subducted SAMW and AAIW decreases due to the buoyancy gain and shoaling of the winter mixed layer. Qu et al. (2020) studied the variability in SAMW using the subduction rate during the Argo period and found that there is a strong correlation between the quasibiennial variability in SAMW and the Southern Annular Mode (SAM). Therefore, the accurate calculation of the subduction rate has become a precondition to investigate climate change from the perspective of ventilation. Through our research, the special circumstances that may cause deviations are taken into account to ensure effectiveness when using the classical method to calculate the subduction rate of these typical water masses. Apart from the AAIW and SAMW, the subtropical water masses mainly refer to STMW. In the South Pacific, the STMW is more specifically known as ESTMW or South Pacific Eastern Subtropical Mode Water (SPESTMW). SPESTMW spreads northwestward from the ventilation region within the subtropical gyre, eventually joining the South Equatorial Current (Wong and Johnson, 2003; Liu and Wu, 2012). SPESTMW is located at 12S°–24°S and coincides with one of the deviation maxima in the South Pacific, so ESTMW is probably overestimated. Due to the limitations of data and the amount of calculation, the situation in the Northern Hemisphere is not taken into account in this research. However, due to the symmetry of the wind field in the Northern and Southern Hemispheres, it can be predicted that vertical pumping still plays a major role in the formation of STMW in the Northern Hemisphere, and the conditions for the two kinds of deviations will also be satisfied in the tropical and subtropical regions of the Northern Hemisphere, so there should be a similar conclusion. In the area where these deviations exist, Cerovečki and Giglio (2016) and Cerovečki et al. (2019) analyzed the T/S using Argo profiles for 2005 to 2012 and related the decrease in volume of STMW to the change in pycnocline in the formation region, which is closely associated with the subduction process. However, the existence of STMW also has a substantial impact on the pycnocline, and the presence of thick STMW will shoal the upper pycnocline and cause subtropical front (STF) variations on decadal timescales in the North Pacific (Kobashi et al., 2021). As a result, because the main formation regions of SAMW and AAIW are at middle latitudes in the Southern Hemisphere, it is accurate to use the classical tracing method to calculate the subduction rate related to these water masses, but the possible deviations should be accounted for when considering subtropical water masses.

As stated above, the Lagrangian method is dealt with directly when we consider the deviation, but the velocity β-spiral is an unavoidable problem. To solve the problem, there is an effective method, that is, to use high temporal resolution data and the Eulerian method. The Eulerian method is accurate in theory and adopted by many previous studies (e.g., De Szoeke, 1980; Cushman-Roisin, 1987; Liu et al., 2016), but it has a large limitation; that is, it requires a very high temporal resolution of data. For most of the subduction areas, the fluctuation of the mixed layer in the selected Stommel's demon period can reach hundreds of meters. Therefore, the deviation caused by the selection of the effective subduction period determined by the monthly averaged data will be very large, so this scheme is not feasible in practice. One possible approach is to use linear interpolation and directly interpolate the monthly averaged data into quasi-daily averaged data. Considering the Lagrangian tracing process in daily steps, the computational cost would be very high, and the implementation of this method needs to be further simplified in follow-up research. In this paper, the Eulerian method is used in the discussion of the second deviation, that is, to judge whether the subduction in the non-Stommel's demon months is effective. The reason we can use the Argo data with low resolution here is that the fluctuation in MLD is no longer intense in these periods, and the amount of subduction in each month is no longer large, so the absolute deviation in judging the effective period is small. At the same time, as a conservative estimate, only the month of fully effective subduction is considered, so the second deviation is the minimum estimation. In particular, it should be pointed out that the correction Equation of comprehensive deviation (Eq.9) is exactly equivalent to the result of using the complete Eulerian method if the location has both the first deviation and the second deviation. Because at this point, the effective subduction time is longer than the selected Stommel's demon period, so the deviation caused by the time resolution is relatively small.

Further analysis shows that the conditions for the two deviations are similar: strong vertical pumping in the source area and along the path of subduction, as well as MLD with little change. Therefore, the distribution areas of the two deviations are similar but not coincident. Due to the effect of vertical pumping, for the first deviation, vertical pumping along the path is the main influencing factor, which forces the water mass subducted at the source to sink to a very deep depth, thus exceeding the amount of water that the source can produce in winter. For the second deviation, both the vertical pumping along the path and local vertical pumping are very important. The former increases the effective subduction time window, while the latter determines the amount of subducted water that can contribute to the effective subduction window. The small change in the mixed layer also reduces the effect of the lateral induction term and makes the effect of the vertical pumping term more significant.

5 CONCLUSION

Using low-resolution monthly averaged observation data, at least two possible deviations are discussed in the classical method: the first is that when the mixed layer becomes deeper during the period of the first deep winter, less water may be subducted at the source, which leads to the overestimation of the subduction rate; the other is that the subduction rate may be underestimated without considering subduction during months other than deep winter months at the source. The results of quantitative analysis using Argo and NCEP data show that the two kinds of deviations exist significantly in the low- and middle-latitude subduction areas, their distribution patterns are similar, and they can partially offset each other, but the comprehensive deviation is still significant and cannot be ignored. At middle latitudes, there are almost no such two deviations. The evaluation of these two deviations has distinct physical meanings, and the application of the revised subduction rate in the global ocean will be a meaningful topic for future study.

6 DATA AVAILABILITY STATEMENT

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGMENT

The authors thank all the reviewers for their comments to improve the original manuscript.