## Abstract

The zonally integrated flow in a basin can be separated into the divergent/nondivergent parts, and a uniquely defined meridional overturning circulation (MOC) can be calculated. For a basin with significant volume exchange at zonal open boundaries, this method is competent in removing the components associated with the nonzero source terms due to zonal transports at open boundaries. This method was applied to the zonally integrated flow in the Indian Ocean basin extended all the way to the Antarctic by virtue of the ECCO dataset. The contributions due to two major zonal flow systems at open boundaries, the Indonesian Throughflow (ITF) and the Antarctic Circumpolar Current (ACC), were well separated from the rotational flow component, and a nondivergent overturning circulation pattern was identified. Comparisons with previous studies on the MOC of the Indian Ocean in different seasons showed overall consistency but with refinements in details to the south of the entry of the ITF, reflecting the influence of ITF on the MOC pattern in the domain. Other options of decomposition are also examined.

## 1. Introduction

The oceanic general circulation is a complex dynamic system in the three-dimensional space. To gain physical insight, some two-dimensional-section views have been widely used. For example, the meridional overturning circulation (MOC) in the ocean has been explored for a long history (e.g., Siedler et al. 2013); in particular, the corresponding streamfunction has been used as a tool to visualize and quantify the MOC.

The common practice in defining the MOC streamfunction is as follows. First, the three-dimensional velocity is zonally integrated to obtain the zonally integrated velocity field (*V*, *W*). Second, the streamfunction is calculated by vertically integrating the zonally integrated meridional velocity. This method has been successfully applied to the world oceans with a particular focus in the Atlantic Ocean (Cabanes et al. 2008; Danabasoglu et al. 2014; Köhl and Stammer 2008; Lumpkin and Speer 2007).

For the Atlantic Ocean (north of 34°S and with the Mediterranean and the Gulf of Mexico included) there is no open zonal boundary; thus, if we omit the freshwater fluxes associated with precipitation/evaporation and river runoff, the zonally integrated flow is approximately nondivergent. But for the general cases, there may be volume flux through the zonal boundaries and the zonally integrated (*V*, *W*) field is no longer nondivergent; thus, the simple method mentioned above is no longer valid.

In fact, there is a large volumetric exchange, the Indonesian Throughflow (ITF), between the Pacific and Indian Oceans. Thus, the zonally integrated velocity field in both the Pacific and Indian Ocean is divergent. Based on the Helmholtz’s decomposition theorem, any two-dimensional velocity vector **V** can be decomposed into an irrotational (divergent or potential) part **V**_{Φ} and a nondivergent (rotational) part **V**_{Ψ}, corresponding to a potential function and a streamfunction, respectively. A straightforward vertical integration of the nondivergent portion **V**_{Ψ} then leads to the MOC streamfunction; with adequate boundary constraints, this streamfunction can be defined unambiguously, and independent on whether the integration is started from the upper or lower boundaries. If the overturning circulation streamfunction is directly calculated by simple vertical integration from the original flow field without such decomposition, the result will be sensitive to the choice of the vertical boundary where the integration starts.

Derivation of the overturning circulation in many parts of the world oceans has similar problems. For example, the zonally integrated velocity field in the South Pacific Ocean has a nonzero source term due to the ITF; the South China Sea has nonzero source term because of the South China Sea Throughflow. More importantly, if one wants to study the meridional circulation in each subbasin in the Southern Ocean, the source term associated with the Antarctic Circumpolar Current (ACC) must be handled carefully as well. With a similar method, Watterson (2001) obtained a MOC pattern of the Arctic and Atlantic which extends all the way to the Antarctic by removing the source term due to ACC at the southern part of the domain. We develop a different simple numerical method to solve the problem in the Indian Ocean where the nonzero source terms involve not only the ACC but also the ITF based on the more up to date dataset.

In this study, we formulate the velocity decomposition problem in rigorous terms mathematically. The fundamental aspects of decomposition, including discussions on the nonuniqueness of the solution and the numerical scheme, are outlined in section 2. The implementation of the method to the Indian Ocean extended to the Antarctic with the dynamically consistent ocean state estimate, Estimating the Circulation and Climate of the Ocean (ECCO), is discussed in section 3. The derived MOC patterns and comparison with previous research are presented in section 4, as well as the discussion on an alternative option of performing the decomposition. Conclusions of this paper are summarized in section 5.

## 2. Method

### a. Decomposition of the zonally integrated flow vector

For an ocean model based on volumetric conservation, the continuity equation is

where (*u*, *υ*, *w*) denote the three-dimensional velocity components. Zonally integrating Eq. (1) in a basin leads to

where *x*_{w} = *x*_{w}(*y*) and *x*_{e} = *x*_{e}(*y*) are the zonal boundaries; *u*_{w} and *u*_{e} denote the zonal velocity at the western/eastern boundaries of the model domain. Denote the zonally integrated flow components as

thus, Eq. (2) can be rewritten as

where *S*(*y*, *z*) is the divergence term, which can be denoted as *S*_{υw} or *S*_{u}. As shown in the lhs and rhs of Eq. (2),

A common practice is to assume that the MOC streamfunction *ψ*_{MOC}(*y*, *z*) satisfies the following relation:

Accordingly, the MOC streamfunction can be calculated by vertically integrating the zonally integrated meridional flow component *V*(*y*, *z*). However, if the source term in Eq. (4) is nonzero, the streamfunction obtained by such a simple vertical integration of the zonally integrated meridional velocity is not uniquely defined. In fact, if we start the integration from the upper or lower boundaries, streamfunction obtained is not zero at the end of the vertical integration. As a result, the corresponding MOC cannot be uniquely defined.

One example of such kind of problem is related to defining the MOC in the Indian Ocean basin. Due to the inflow of the ITF at the eastern boundary of the basin, there is a net meridional volume transports at all latitudes south of the ITF (starting from around 10°S). Since the meridional transport associated with the ITF is mostly confined to the upper part of the water column, in the study of the deep cell associated with the MOC, investigators chose to start the integration from the lower boundary, trying to minimize the impact of the ITF on the deep MOC cell (Wang et al. 2014, hereafter WZWK). However, the impact of the ITF is not entirely confined to the upper ocean as will be seen in the following discussion; thus, such approach may not provide an accurate description of the MOC cells. In this paper, we aim to provide a definition of the overturning circulation streamfunction that is less dependent on the choice of boundary from which the integration Eq. (7) is started.

In the general case, the velocity vector of the ocean circulation can be separated into an irrotational (divergent or potential) part and a nondivergent (rotational) part according to the Helmholtz decomposition (e.g., Arfken and Weber 2005). The Helmholtz decomposition theorem has been applied in electrodynamics at first (e.g., Griffiths 2013), and then also widely adopted in meteorology in decomposing the wind field into rotational and divergent parts for the atmospheric circulation analysis (e.g., Hawkins and Rosenthal 1965; Sangster 1987; Knutson and Weickmann 1987; Frederiksen and Frederiksen 1993) and data assimilation research (e.g., Barker et al. 2004; Rawlins et al. 2007). The decomposition method has also been applied in the ocean circulation problems in decomposing the surface currents (Li et al. 2012), horizontal gyre, the overturning circulation (Watterson 2001), and ocean heat transport in ECCO (Forget and Ferreira 2019).

In our case, the zonally integrated, two-dimensional velocity field **V** with components as in Eq. (3) is the targeted vector and can be decomposed as (e.g., Watterson 2001; Arfken and Weber 2005; Li et al. 2012)

where **V**_{Φ} and **V**_{Ψ} represent the irrotational (or divergent) part and the nondivergent (or rotational) part with the following definitions, respectively,

where scalars Φ and Ψ denote the potential function of the irrotational component and the streamfunction of the nondivergent component of the original flow vector **V**, and **k** denotes the unit vector in the (eastward) zonal direction.

In the scalar form, the decomposition of the 2D flow field (*V*, *W*) can be written as

where

in which *V*_{Φ} and *W*_{Φ} refer to the irrotational components, and *V*_{Ψ} and *W*_{Ψ} the nondivergent components of the flow field (*V*, *W*) respectively.

This is a 2D Poisson equation for the potential function of the irrotational component. The abovementioned decomposition can separate the irrotational component **V**_{Φ} from the total meridional transport, and this component corresponds to the net zonal inflow at the open boundaries, Eq. (6). By subtracting **V**_{Φ} from the total flow vector **V**, one obtains a clean nondivergent **V**_{Ψ}, which provides a unambiguous overturning circulation.

### b. The nonuniqueness of the decomposition and boundary conditions

The Helmholtz decomposition into the rotational and divergent components is unique for infinite domain. However, for a limited domain, the decomposition is not unique due to lack of additional physical constraints on the boundary conditions (BCs) for either of the decomposed component. This point has been discussed in both meteorology (Bijlsma et al. 1986) and oceanography (Fox-Kemper et al. 2003). The global domain with periodic domain is free of such problem but additional assumptions on boundary conditions have to be made if one wishes to obtain a solution of decomposition for a problem with presence of boundary of the computational domain. Many attempts have been made in meteorology [see Li et al. (2012) for a list of those studies] and oceanography (Bryan et al. 1999; Roberts and Marshall 2000; Watterson 2001) to find the appropriate additional constraints; these solutions represent choices of *a* decomposition, not *the* decomposition, as elaborated in Fox-Kemper et al. (2003).

We acknowledge the nonuniqueness of the Helmholtz decomposition in a limited domain, and will work on obtaining a decomposition with the boundary conditions consistent with the aim of this paper, that is, to derive an overturning circulation streamfunction that is independent on the boundaries from which the integration Eq. (7) starts. This requirement leads to the first constraint, that is, the net meridional transport of the nondivergent component across any latitudinal section vanishes or

where *L* denotes the length of the domain in *y* or meridional direction, and *H* the depth of the ocean floor, which is allowed to be function of *y*. And the second constraint is the same as being adopted by Watterson (2001), that is, the streamfunction Ψ is constant along wall boundaries, or in another word, the curve along the wall behaves as a streamline,

The wall boundaries include the ocean floor at bottom and the sidewalls at both ends of the domain in the meridional direction. Note that the wall boundaries come from zonal integration of the ocean domain, thus depend on the arbitrary choice of the domain of computation (see Fig. 1 for an illustration).

Equations (15) and (16) comprise all the assumptions we shall make throughout this study. Boundary conditions for the potential function Φ and the streamfunction Ψ will all be derived in the following from these two assumptions.

With Eqs. (13), (15), and (16), the surface boundary condition for the streamfunction Ψ can be readily obtained as

which implies that the free surface behaves as a streamline with the same value as that along the wall boundary. In another word, the curve of the constant streamline encloses the whole computational domain in the latitude–depth plane.

Next will be the boundary conditions for the potential function Φ. With the nonpenetrating condition for the original flow field **V** and Eqs. (8)–(10), we can obtain

where **n** denotes the normal unit vector of the wall boundary. Since constant streamfunction is assumed, the second term in above equation vanishes due to

with **l** the tangential unit vector of the wall boundary, which yields the bottom boundary condition for Φ as

At surface, since Ψ keeps constant, the vertical motion of the rotational part *W*_{Ψ} vanishes according to Eq. (13). Thus as per the velocity decomposition in Eq. (11), the surface boundary condition for Φ, takes the form

If the model is based on volume conservation (such as the ECCO product to be introduced in the next section), *W*|_{z=0} in above equation represents the net volume flux across the sea surface (evaporation minus precipitation). Therefore, under the two assumptions made previously, the boundary conditions for both Φ and Ψ are complete. The system of Poisson equation for the potential function Φ can be formed as below:

where Ω denotes the latitude–depth computational domain. Watterson (2001) proposed an efficient iterative method of eliminating the spurious vorticity at boundary to solve the Poisson problem of the streamfunction Ψ in the limited domain with a standard code of Poisson solver. The MOC streamfunction of the Atlantic (including its Southern Ocean sector) and Arctic were successfully obtained with his method from a model’s data. In this paper, we use a different method of incorporating the boundary conditions (both surface and bottom) into the numerical iteration of the interior points, and solving the finite-differenced equation of the potential function Φ as a whole. The nonuniform grids of the ECCO data are appropriately transformed to be uniform with some techniques, taking advantage of the best performance of the numerical scheme under uniform grids. This method is straightforward in physical meaning and simple to implement. And the details are explained in the following section.

On the other hand, for a Poisson equation with purely Neumann BCs as in Eq. (21), there is no solution unless a compatibility condition is satisfied (e.g., Li et al. 2017), which is

where ∂Ω denotes the boundary of domain Ω, *dA* and *ds* the areal element and line element of domain Ω and boundary ∂Ω, respectively, and (∂Φ/∂*n*) the outward normal derivative of Φ along boundary. With the boundary conditions of Φ as shown in Eqs. (21) and (22) can be rewritten as

This compatibility condition will be checked in section 4 for the reliability of the numerical solution.

It is also worth mentioning that if the compatibility condition is satisfied and ∂Ω is smooth, the solution does exist but it is not unique. Indeed, for any solution Φ_{0}, Φ_{0} + *C*_{0} (*C*_{0} is an arbitrary constant) is also a solution (e.g., Li et al. 2017). However, our concern in the flow decomposition problem is limited to the velocity component (*V*_{Φ}, *W*_{Φ}), which are unique, regardless any additive constant for Φ.

### c. Numerical solution to the Poisson equation

In general, the two-dimensional Poisson equation [Eq. (21)] can be solved numerically using the “successive over-relaxation iteration” (SOR) method, which converges faster than either Jacobi or Gauss–Seidel methods (e.g., Li et al. 2017). For a uniform 2D mesh, the general scheme of the SOR method is as follows

where the relaxation parameter *ω* determines the iterations needed for convergence in different setups, and is usually chosen within the range 1 ≤ *ω* < 2. The superscripts *n* − 1 and *n* denote the previous and the current iterative step, respectively. The subscripts *j* and *k* denote the grid point index in *y* and *z* directions of the discretized mesh, respectively. The term Δ*h* denotes the uniform grid size in either direction.

With regards to the incorporation of the BCs, it can be readily seen that Eq. (24) can be modified as the following form with Neumann BCs at all boundaries,

where $Aj,k$ denotes the mask matrix recognizing the numerical grid point as water or nonwater,

and

represents the total number of water grids adjacent to grid (*j*, *k*); the $Tj,kn$ term in Eq. (25) is introduced in order to deal with the nonhomogeneous Neumann BC at the surface,

where *W*_{j,1} = (Φ_{j,1} − Φ_{j,2})/Δ*h* = *W*_{surface} is the vertical derivative adjacent to the surface; this surface vertical velocity represents the volume flux across the sea surface. Abovementioned numerical schemes are all constructed based on the Arawaka C grid in consistency with ECCO.

In this study our testing indicated that *ω* = 1.995 gives the fastest convergence. The optimal value of *ω* adopted here is in good agreement with the result of the classical theory on the iterative methods (Yang and Gobbert 2009). Other options of *ω* can also produce the same iterative results, but with a slower convergence of the runs.

Therefore, Eqs. (25)–(28) comprise the numerical solution system to the 2D Poisson problem, Eq. (21). The surface BC and the irregular bottom topography is implicitly implemented through the specification of the mask matrix. In the following, this method will be applied in decomposing the flow components for the sake of deriving a consistent MOC pattern in the Indian Ocean extended to the Antarctic.

## 3. Data and implementation

As discussed above, the zonally integrated velocity field in the Indian Ocean basin extended to the Antarctic has nonzero source terms due to the ITF and the ACC, which impedes derivation of a physically consistent MOC pattern directly from original flow field. Thus, we apply the method developed in section 2 to separate those source terms from the original flow field before deriving the overturning streamfunction.

The only input variables we need to achieve this goal is the 3D flow field in the Indian Ocean. The dataset we rely on is the ECCO state estimate release 4, version 3 (v4r3); this dataset can be interpreted as the result of the least squares fit of the MITgcm (Adcroft et al. 2004), with sea ice and mixed layer submodels to about *O*(10^{9}) observational points (Forget et al. 2015; Fukumori et al. 2017). Lagrange multipliers are used (Thacker and Long 1988) to enforce the model equations in such a way that basic conservation rules (heat, freshwater, momentum, energy, etc.) are satisfied globally and locally for the period of January 1992–December 2015. In the horizontal directions, ECCO v4r3 adopts a nonuniform LLC90 (latitude–longitude–cap) gridding system, which maps the earth using five faces including an Arctic face and four mostly latitude–longitude sectors, with a nominal resolution of 1°. In the vertical direction, ECCO has rescaled height coordinates with 50 vertical levels and partial cell representation of bottom topography (Campin et al. 2008). More information on ECCO is available in Forget et al. (2015) and on the website ecco.jpl.nasa.gov.

In extracting the flow field for the Indian Ocean, the ECCO dataset provides the basin mask. To accommodate with our application, we slightly modify its open boundary at the ITF by moving it 8° west in order to turn the meridional throughflow via Lombok strait to be incorporated into the zonal mainstream of the ITF (Fig. 1a). With this treatment, all the sources of flow divergence as in Eq. (4) reside in the zonal inflow/outflow at the open boundary. The domain we choose for computing the streamfunction of the MOC in the Indian Ocean is illustrated in Fig. 1. This domain spans 20°–147°E in longitude and 70°S–26°N in latitude, and it involves three faces (faces 1, 2, and 4) in the LLC gridding system introduced in ECCO v4 (Forget et al. 2015).

In previous studies, the southern boundary of the Indian Ocean is often posed around 34°S to exclude the impact of the strong zonal flow, the ACC (e.g., Lee and Marotzke 1998, hereafter LM; WZWK). However, such a southern boundary requires additional information about the irrotational flow component **V**_{Φ} at this open meridional boundary when solving the Poisson equation, Eq. (21). It is clear that such a boundary condition is not available. To overcome this problem, we set the model domain as the whole Indian Ocean, including the sector of the Southern Ocean all the way south to the Antarctic Continent. Regarding the time resolution of the ECCO product, we adopt the monthly mean data, which include a total of 288 snapshots from 1992 to 2015 for each time-dependent variable.

The grid size of LLC in meridional direction Δ*y* is not only nonuniform in latitude, but also varies from longitude to longitude, which is illustrated in Fig. 1c. However, within the latitudinal range of our computational domain, that is, 70°S–26°N (Fig. 1b), Δ*y* does not change with longitude (Fig. 1c). This facilitates the operation of zonally integrating the flow field along certain longitude, which yields (*V*, *W*) with Eq. (3). This is the targeted flow vector **V** that we shall decompose with the method introduced in section 2. It is worth mentioning that the flow vector adopted in this paper refers to the Eulerian velocity; thus, the bolus velocity which is a parameterization of the mesoscale eddies (Gent and McWilliams 1990) is not included in our MOC computation though it is available in ECCO. As a result, the MOC derived in this paper is a Eulerian one. Regarding the surface boundary condition in Eq. (20), the freshwater flux across the sea surface (evaporation minus precipitation) available as a standard ECCO diagnostic is applied in our calculations.

Before solving the Poisson equation [Eq. (21)], the source term *S* should be calculated. As an example, we randomly pick one month (January 1997) and plot the corresponding terms associated with the divergences of *V* and *W*. Figure 2 shows that *S*_{υ} and *S*_{w} cancel each other very well at most latitudes, except in the latitude bands of the open zonal boundaries, ITF (around 20°–10°S as can be seen from Fig. 1a) and ACC (south of 34°S), which simply indicates the nonzero volume exchange across the open boundary. The other conclusion that can be drawn from Fig. 2 is that the feature of *S*_{υw} in ACC below 100 m is mainly attributed to *S*_{υ}.

Figure 3 compares *S*_{υw} and *S*_{u}, which shows that these two terms are almost equivalent to each other over the whole domain, with a root-mean-square (RMS) as small as *O*(10^{−4}) m s^{−1}, two orders of magnitude smaller than *S*_{υw} or *S*_{u} per se. The nearly identical outcomes produced by two independent methods verify a reliable source term and reflect the outstanding conservation property of the ECCO product as well to a large extent.

Another technical issue with the latitude/depth grids of ECCO is that it is nonuniform in both directions while all the tests and discussions in section 2 are performed with uniform grids. Even though a SOR scheme for nonuniform grids can be developed, the problem associated with the huge width-to-height aspect ratio, that is, Δ*y*/Δ*z* = *O*(10^{3})–*O*(10^{4}) in the original ECCO data, makes the iterative process infinitely long and the iterative performance substantially degraded in practical application. This issue is confirmed in an idealized numerical test (not shown here). To overcome this problem, we propose a workaround with the following steps:

Interpolate the source term from the original nonuniform ECCO grid to a new set of grid that is uniform in each direction, 10 km in the

*y*direction and 10 m in the*z*direction for example. The interpolation method must be integral preserving so that the domain integral of the source term is kept unchanged before and after the interpolation.Rescale the meridional grid with a factor of

*δ*= 1/1000. In this way, both the meridional coordinate and the meridional velocity is rescaled with the same factor of*δ*; in this new meridional coordinate, $y*=\delta y$ and $V*=\delta V$; thus this rescale leaves the divergence contribution due to the meridional velocity*S*_{υ}unchanged.

The model domain is now covered by a new computational mesh with square grids of 10 m × 10 m. Thus the Poisson equation for the potential function $\Phi *$ under the transformed coordinate satisfies

To summarize the procedure above, the nonsquare mesh problem can be overcome by transforming to the square mesh ($\Delta y*=\Delta z$). The new set of the Poisson equation bears the same form as the original one, Eq. (21), except with Φ and *y* being replaced by $\Phi *$ and $y*$ and that the way of recovering the meridional flow component now becomes

With regards to the computational domain for our Indian Ocean case, as being pointed out earlier in this section, both meridional borders are placed in land with $y*=0$ corresponding to 70°S and $y*=L$ corresponding to 26°N; the choice of the vertical boundaries are rather straightforward with *z* = 0 at the surface and *z* = −*H*(*y*) at the ocean floor.

Now that the nonuniform grids can be transformed as square with above steps, the SOR scheme for square grid, Eq. (24), can be applied directly to the ECCO data with no necessity of developing a nonuniform version of the iterative scheme.

With the procedures discussed above, the set of equations, Eq. (21), can be solved numerically by the modified SOR iterative method for a limited domain with the ECCO data to yield an irrotational component (*V*_{Φ}, *W*_{Φ}) in the domain of the Indian Ocean extended to the Antarctic indicated by Fig. 1a. The stopping criterion for the numerical iterative process chosen for this case is the RMS difference of $\Phi *$ between two consecutive iterative steps is less than 10^{−8}. The solution and further results on the cells of MOC in the Indian Ocean extended to the Antarctic are presented and discussed in the next section.

## 4. Results

### a. The Eulerian MOC of the Indian Ocean

Immediately after the irrotational component (*V*_{Φ}, *W*_{Φ}) is derived by solving the Poisson problem in section 3, the nondivergent component (*V*_{Ψ}, *W*_{Ψ}) can be derived by subtracting it from the total flow (*V*, *W*) as per the decomposition in Eq. (11). As assumed in Eq. (15), the vertical integration of *V*_{Ψ}(*y*, *z*), that is, the net meridional transport by the nondivergent flow component, must vanish at all latitudes. However, the numerical algorithms in our flow-decomposition method unavoidably introduce errors which make the net meridional transport by *V*_{Ψ}(*y*, *z*) deviate from zero. In this sense, the net meridional transport by *V*_{Ψ}(*y*, *z*),

may serve as a measure of the decomposition error. The net meridional volume transports for the total flow *V*, the irrotational component *V*_{Φ}, and the nondivergent component *V*_{Ψ} are plotted in Fig. 4 for January 1997 (Fig. 4a) and July 2010 (Fig. 4b), two randomly picked snapshots from the ECCO time series for demonstration. One can find that the net transports by *V* and *V*_{Φ} are very close to each other, and the corresponding error $M\Psi net\u2061(y)$ is two orders of magnitudes smaller than the magnitude of net transports by *V* and *V*_{Φ} at all latitudes. The RMS values of $M\Psi net\u2061(y)$ for the two cases in Fig. 4 are 0.100 and 0.113 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}), respectively. The low level of the decomposition error justifies the effectiveness of our separating algorithm.

Another way to assess the effectiveness of the decomposition is through the divergence and curl fields of the decomposed nondivergent and irrotational components. Using a randomly selected month (January 1997) from ECCO as an example, the divergences and curls of the three flow fields are plotted in Fig. 5, respectively. It is evident that virtually all the divergences of the original flow field are taken by its divergent or irrotational part, leaving the nondivergent or rotational part well divergence-free; and the curls of the original flow field are taken almost all by its nondivergent or rotational part, leaving the divergent or irrotational part well curl-free. Using the root-mean-square value as the measure, the divergent part and the rotational part account for around 98.8% of the divergence and 99.8% of the curl of the original flow field, respectively. Thus, our decomposition method has done a fair job in obtaining the nondivergent and irrotational components of the original flow field.

The flow decomposition is to obtain a physically consistent pattern of the Eulerian MOC, which is the aim of this paper. As a common practice, the streamfunction of the MOC is derived by vertically integrating the meridional flow component accumulatively in the vertical direction either upward or downward by assuming zero streamfunction at starting point, with the following forms consistent with Eq. (7):

where $\psi MOCup$ and $\psi MOCdown$ refer to the upward and downward integration, respectively.

To demonstrate the effect of the flow decomposition on the patterns of MOC, we calculate the streamfunctions with both Eqs. (32) and (33) for the total meridional component *V*, its irrotational part *V*_{Φ}, and nondivergent part *V*_{Ψ}, respectively. These streamfunctions are named *ψ*(*V*), *ψ*(*V*_{Φ}), and *ψ*(*V*_{Ψ}) and are demonstrated for January 1997 (Fig. 6) and July 2010 (Fig. 7). One can see that the integration directly with the total meridional component *ψ*(*V*) (Fig. 6 and Fig. 7, upper row), cannot provide a consistent MOC pattern in the domain with MOC cells strongly depending on the direction or starting point of the integration, in both the ACC and the ITF regions where significant zonal inflow/outflow occurs. Nevertheless, if the irrotational part *ψ*(*V*_{Φ}) (Fig. 6 and Fig. 7, middle row) that accounts for the zonal flow impact is removed from the total, yielding the nondivergent part *ψ*(*V*_{Ψ}) (Fig. 6 and Fig. 7, bottom row), the patterns of the Eulerian MOC become consistent between the two ways of integration. In July (Fig. 7), two MOC cells are well defined by the nondivergent streamfunction *ψ*(*V*_{Ψ}): a clockwise cell (positive streamfunction) in the ACC region above 500 m and a large-scale counterclockwise cell (negative streamfunction) prevails in the whole depth of the rest of the basin, with a significant part existing in the Indian Ocean basin (north of 30°S). While in January (Fig. 6), a new clockwise cell is witnessed in the northern area and becomes the dominant MOC cell in the Indian Ocean. The MOC pattern looks rather stable throughout the year in the ACC region, while it exhibits notable seasonal variability in the Indian Ocean.

Though the results of only two snapshots were demonstrated above, they are representative among the total 288 monthly fields of ECCO v4r3 from January 1992 to December 2015. For an overview of the effectiveness of the method throughout the ECCO’s time series, the mean value and standard deviation of the error indicator $M\Psi net\u2061(y)$ for each month are presented in Fig. 8a, with the two vertical dashed lines indicating January 1997 and July 2010, whose MOC patterns have been illustrated previously in Figs. 6 and 7. The errors, mostly within 0.1 Sv with a standard deviation of around 0.1 Sv, show negligibly small values overall, compared with the much larger net transport of the total flow, whose RMS reaches around 30–40 Sv. This uniformly small value of $M\Psi net\u2061(y)$ also proves that the previous assumption, Eq. (15), is kept valid in solving the problem. And the well constrained residual transport warrants the derivation of a physically consistent MOC streamfunction. On the other hand, the compatibility condition of the existence of solution for the Poisson problem mentioned in section 2, Eq. (23), is also checked, with result shown in Fig. 8b. It is apparent that the areal integral of the source term [lhs of Eq. (23)] and the line integral of the zonally integrated vertical velocity at surface [rhs of Eq. (23)] well equate with each other throughout the time series of the ECCO data, with an RMS difference as small as 1.5 × 10^{−4} Sv. The balance of these two terms implies that the rate of change of the oceanic volume in the domain of interest (Fig. 1a) is entirely due to the volume exchange at open boundaries in zonal. And the satisfaction of this existence condition warrants the credibility of the solution derived in this section.

Previous research on the pattern of the MOC cells in the Indian Ocean were mostly confined to the south of ACC, or south of around 34°S to exclude the influence of the strong zonal volume exchange (LM; Ferron and Marotzke 2003, hereafter FM; WZWK). However, the influence of the ITF still remains in the domain. To better compare with those results, we also narrow the region to the same domain and present climatological streamfunctions of the MOC for January (Fig. 9) and July (Fig. 10). The area between the vertical yellow lines indicate the latitudinal range of ITF at the eastern boundary of the basin. In fact, the ITF enters the Indian Ocean from the eastern boundary, with a big portion joining the South Equatorial Current and flowing southward after reaching the western part of the Indian Ocean basin (WZWK). As a result, the pattern of *ψ*(*V*) is very similar to that of *ψ*(*V*_{Ψ}) north of the ITF entry in both months. However, notable differences do emerge within and south of the ITF entry, especially near the surface, manifesting the impact of the ITF and its down stream.

The climatological pattern of the Indian Ocean MOC in January and July is also shown by WZWK (their Fig. 1) and LM (their Fig. 20); their figures exhibit a similar overall pattern and value of the MOC as our results in Figs. 9 and 10. The details of the MOC show better similarity north of the ITF than south since there is less zonal flow impact in the north [the volume transport of the Malacca Strait is only an order of magnitude of 0.1 Sv, according to SODA and ROMS (Daryabor et al. 2015)]. But when it enters the ITF region, some differences do occur. For the January case, the counterclockwise cell in WZWK is limited to a rather shallow depth north of 20°S; the same cell of LM does reach 500 m similar to our result, but it failed to extend downward all the way to the bottom as our and WZWK’s results show. For the July case, our and LM’s results are both dominated by a widespread counterclockwise cell in the basin, while in that of WZWK, a strong clockwise cell prevails in the whole water column between 34° and 20°S. Such a pattern is quite similar to the streamfunction integrated by the total flow component *V* in our Fig. 10a. Actually, the ITF’s influence on the pattern of the MOC was realized by WZWK who tried to minimize it by integrating “the mass streamfunction from the ocean bottom upward such that ITF’s influence is confined within the upper ocean since most of the ITF transport is present in the upper around 400 m.” But even so, the ITF’s influence on the MOC is still remarkable as seen in above comparisons. From Figs. 9b and 10b, one can find the ITF influence as a source term is not limited within the upper layers, but can reach several thousand meters deep, especially in boreal summer. In fact, notable eastward volume transport exists at depths underneath the westward mainstream of the ITF at the eastern entrance of the domain according to ECCO (Figs. 9d and 10d). The eastward transport below 400 m is estimated to be around 15%–32% of the total flux of the westward transport for various seasons.

FM worked out an annual mean MOC in the Indian Ocean from a 4D-variational assimilation of WOCE hydrography in the year 1995 (their Fig. 6b). When it is compared with the decomposed annual mean MOC of the same year by our method based on the ECCO data (Fig. 11), one finds that the overall pattern of the anticlockwise cells is rather similar, with a surface cell of 10 Sv (FM) and 16 Sv (this study) and a deep cell of 18 Sv (FM) and 14 Sv (this study), respectively. However, a clockwise cell above 500 m around 30°S remains in FM. It resembles the pattern of MOC derived without decomposition (Fig. 11b) reflecting the impact of ITF that is supposed to be removed from the MOC. It is also worth noting that the annual mean MOC without decomposition (Fig. 11a) produces too weak deep cell of only 4 Sv. This result throws some doubts on the attempt of isolating the impact of ITF on the deep cells by integrating the MOC streamfunction upward from the bottom (WZWK).

In summary, in our results, a Eulerian MOC can be clearly identified thanks to the rather clean removal of the ITF’s influence with a physically meaningful method.

### b. On the alternative option of flow decomposition

Besides solving the Poisson equation of the potential function Φ, the flow decomposition could also be achieved by solving the Poisson equation of the streamfunction of the nondivergent component. Actually, the system of the Poisson equations for the streamfunction Ψ can be readily formed as below with its BCs set as in the previous section, Eqs. (16) and (17),

Note that the source term of Ψ’s Poisson equation is the curl of the original flow field instead of the divergence as in Φ’s equation. In the application of meteorology, it was stated that the decompositions by solving Φ or Ψ are equivalent (Li et al. 2012). To further verify the results from our method, we have explored the possibility of performing the decomposition by virtue of Ψ’s equation as well, and obtained a clean enough decomposition of two parts, with virtually 100% of the divergence taken by the divergent or irrotational part and over 99.9% of the curl taken by the nondivergent or rotational part. However, the Ψ method is found producing artificial vertical flow components in *W*_{Φ} and *W*_{Ψ}, and the magnitudes of both terms are even notably larger than the original vertical component, *W* itself. Such large errors are related with the large aspect ratio of the computational grid of SOR iterations because the coordinate stretching transformation is not applicable to the Ψ method. Due to the length limit, the details of this problem will not be discussed in this paper.

## 5. Summary and conclusions

The source terms due to exchange through the open zonal boundary impedes the derivation of a physically consistent overturning circulation pattern from the original flow field directly. In an attempt to solve this problem and obtain the streamfunction of the overturning circulation free of such influence, we decompose the laterally (zonally for the meridional overturning circulation) integrated flow vector **V** to an irrotational (divergent or potential) vector field **V**_{Φ} and a nondivergent (sinusoidal or rotational) vector field **V**_{Ψ} as per the Helmholtz’s decomposition theorem. It is known the Helmholtz decomposition is not unique for a limited domain, and additional constraints must be imposed on the boundary to select a solution of decomposition. Compatible with the goal of this paper, two reasonable constraints have been proposed to provide boundary conditions at ocean floor and free surface. The first one is constant streamfunction along the wall boundary, and the second one is the vanished net meridional transport across any vertical section for the rotational (nondivergent) flow field. The first constraint is a commonly adopted one, while the second is necessary for the existing of an overturning circulation streamfunction with no dependence on the starting point of the vertical integration.

A modified successive overrelaxation (SOR) scheme has been developed to solve the set of the Poisson equation for the potential function in the limited domain, with the competence in resolving the boundary conditions at both the surface and the ocean floor with irregular topography. The numerical algorithm we developed is straightforward and simple to implement.

This modified SOR method was applied to the Indian Ocean, which is extended to the Antarctic. The ECCO product (v4r3), a dynamically consistent ocean state estimate dataset based on MITgcm, has been used for this study because of its outstanding conservation property.

The nondivergent flow field **V**_{Ψ} was derived by solving the Poisson equation of the potential function. It was demonstrated that the streamfunctions integrated with **V**_{Ψ} could provide an unambiguous pattern of the Eulerian MOC cells in the domain of study under the settings of the problem, reflecting very little dependence on the direction/beginning of the integration in the vertical. These satisfactory results stem from the near-perfect cancelling of the net meridional volume transports between **V** and **V**_{Φ}, which gives rise to close-to-zero net transports by **V**_{Ψ} at all latitudes, and the fairly clean separation of the divergence and curl of the original flow into the irrotational and nondivergent parts, respectively. The alternative way of performing the decomposition by directly solving the Poisson equation of the streamfunction Ψ was also investigated. It was found that such an approach can lead to artificial vertical velocity of large amplitude in both of the decomposed parts in the interior of the domain, and such artificial velocity would produce a misleading MOC pattern. As a result, although the Ψ method and the Φ method might be equivalent in decomposition for square or near-square grids (such as in the case of the horizontal flow field), they are not for data with a large aspect ratio (such as in the MOC case of this paper).

Compared with that derived by the conventional way of directly integrating the total meridional flow **V**, the Eulerian MOC derived by **V**_{Ψ} sees notable improvements, especially in zones with strong normal flows at zonal boundary. The MOC in the Indian Ocean basin exhibits significant seasonal variability, qualitatively consistent with previous researches but different in details, especially in the zones downstream the ITF. The results also show that the pattern of the deep cells in the Indian Ocean might still be influenced by the ITF even though this inflow is concentrated to the upper levels when entering from the east.

The error of the flow decomposition method, measured by the net meridional volume transport of the nondivergent component, was shown to be negligibly small overall, which warrants the derivation of a physically consistent MOC streamfunction.

We conclude that the method developed in this paper can provide a physically meaningful streamfunction and more consistent pattern of the overturning circulation than the direct integration of the total meridional flow, especially for a domain with sizeable inflow at the lateral boundaries. This decomposition method is compatible with realistic bottom topography. However, there is also limitation with this method: the domain must have a closed boundary (wall or land) at both ends of the meridional direction because neither of the decomposed components is available at the open boundary.

With the approach developed in this paper, the characteristics of the MOC cells, such as their seasonal, interannual, and longer timescale variability, can be investigated. This approach can be applied in studying the overturning circulation in other basins in future. It is also worth mentioning that this method is also applicable to the zonal overturning circulation (ZOC) discussed by Watterson (2001) and Tan et al. (2015).

## Acknowledgments

LH was supported by the National Basic Research Program of China through Grant 2019YFA0606703 and “The Fundamental Research Funds of Shandong University” (2019GN051). The authors thank the anonymous reviewers and the editor for their constructive comments. Code availability: The Matlab code that performs the decomposition and produces some figures in this paper is available at https://github.com/lei-han-SDU/IMOC/.

## REFERENCES

*Mathematical Methods for Physicists*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Ocean Modell.*

*Ocean Modell.*

*Ocean Sci. Discuss.*

*Deep-Sea Res. II*

*Nat. Geosci.*

*Geosci. Model Dev.*

*J. Phys. Oceanogr.*

*J. Atmos. Sci.*

*J. Phys. Oceanogr.*

*Introduction to Electrodynamics*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Appl. Math. Mech.*

*Numerical Solution of Differential Equations: Introduction to Finite Difference and Finite Element Methods*

*J. Phys. Oceanogr.*

*Quart. J. Roy. Meteor. Soc.*

*J. Geophys. Res.*

*Mon. Wea. Rev.*

*Ocean Circulation and Climate: A 21st Century Perspective*. International Geophysics Series, Vol. 103, Academic Press, 904 pp

*Acta Oceanol. Sin.*

*J. Geophys. Res.*

*J. Climate*

*J. Atmos. Oceanic Technol.*

*Appl. Math. Lett.*

Proc. ECMWF Seminar Series on Numerical Methods: Recent Developments in Numerical Methods for Atmosphere and Ocean Modelling, Reading, United Kingdom, European Centre for Medium-Range Weather Forecasts, 139–149, https://www.ecmwf.int/node/7642