Simple algorithms to compute meridional overturning and barotropic streamfunction on unstructured meshes

Computation of barotropic and meridional overturning streamfunctions for models formulated on unstructured meshes is commonly preceded by interpolation to a regular mesh. This operation destroys the original conservation, : which can be then artificially imposed to make the computation possible. An elementary method is proposed that avoids interpolation and preserves conservation in a strict model sense. The method is described as applied to the discretization of the Finite volumE 25 Sea ice – Ocean Model (FESOM2) on triangular meshes. It however is generalizable to collocated vertex based discretization on triangular meshes and to both triangular and hexagonal C-grid discretizations.


Introduction
Over recent years, a considerable progress has been achieved in the development of global ocean circulation models working on horizontally unstructured meshes such as FESOM1.4, Wang et al. (2014), MPAS-Ocean Ringler et al. (2013), FESOM2, Danilov et al. (2017) and ICON-Ocean Korn (2017). By refining in dedicated areas of the world ocean these models may 10 resolve dynamics that would otherwise require nesting or using higher resolution globally. Since these models still use vertically aligned meshes the overhead of horizontally unstructured mesh is minimized because the horizontal neighborhood information is valid for the entire vertical column, and becomes negligible as the number of vertical levels is increased. These models show a very good parallel scalability and reach throughput (in simulated years per day) comparable to that of structured-mesh All new large-scale ocean models are based on the finite-volume method and as such have a clear definition of fluxes at boundaries of the control cells of their meshes. However, these fluxes are defined on irregularly located faces, so instead of using them in their original sense one is tempted to rely on interpolation to a regular mesh. Our practice shows that incurring inconsistencies can be large, and this road should not be followed if global or basin-scale quantities are computed. It turns out that there are efficient and easy-to-implement procedures that are based on exact fluxes and balances and might be used 5 for analyses. These procedures do not rely on interpolation, but use binning, which is sufficient in most cases except for very coarse meshes. The intention of this note is to describe some of them. In doing so we will use the arrangement of variables of FESOM2, however the adjustments needed for other models with different discretizations are relatively straightforward and will be briefly mentioned. We suspect that similar procedures are already used by other groups (in particular, for the analysis on cubed-sphere meshes of the Massachusetts Institute of Technology general circulation model see eg. Adcroft et al. (2004), 10 tri-pole grid configurations of Parallel Ocean Program (POP; Smith and Gent (2004)) or MPAS-Ocean), but we feel that they need to be documented for unstructured meshes, facilitating the use of unstructured-mesh models by a broader community.
We will discuss computations of meridional overturning streamfunction in height and density coordinates as well as computations of barotropic streamfunction.
2 Geometry of discretization 15 FESOM2 uses a cell-vertex discretization, placing horizontal velocities on centroids of triangles and scalar quantities at vertices if viewed from the surface, as shown schematically in Fig. 1. These quantities are stored at midlevels. Vertical velocities are located at vertices and full levels. We use index v to enumerate vertices, c (cells) to enumerate triangles and k to enumerate vertical levels or midlevels (centers of layers). The velocity control volumes are mesh triangles, and scalars are associated with median-dual control volumes formed in the horizontal plane by connecting midpoints of edges with cell centroids. On uniform 20 equilateral meshes they coincide with hexagons of the dual mesh, but in general case they differ. For the reasons discussed in Danilov et al. (2017) the bottom topography of FESOM is given on cells, implying that velocity control volumes are triangular prisms in 3D. However, a part of scalar control volume can be cut by bottom topography at depths, and its footprint will differ from that at the surface. As a consequence, there is a one-dimensional array A c of triangle areas, and a two-dimensional array A kv of the areas of scalar control volumes. The transport through the top face of scalar prism with indices (k, v) is w kv A kv , 25 with w kv the respective vertical (or cross-level in the case of moving level surfaces) velocity. Each triangle is characterized by Fig. 1.
The elementary structure used in computations of horizontal fluxes between two scalar control volumes is given by mesh edges (labelled with index e). An edge is characterized by its two vertices (v 1 , v 2 ) symbolically written as v(e), and two cells it belongs to, (c 1 , c 2 ) symbolically written as c(e). For boundary edges c 2 is absent, and c 1 is the left cell to the edge direction, where e z is a unit vertical vector, h kc1 and h kc2 are the layer thicknesses at respective velocity points and T e is the tracer estimate at edge midpoint. T e = 1 for volume transport. In MPAS-Ocean or ICON-Ocean codes, which are based on hexagonal and triangular C-grid discretizations, normal velocities are located at edges and computations of transports are simpler. The ar-5 rangement of hexagonal C-grid is easily obtained from the case considered here if edges of dual triangular mesh are considered (with the difference that centroids are replaced by circumcenters and lines connecting c 1 with c 2 are perpendicular to edge e).
Importantly, edge-related transports are the same as in the model, however care should be taken that also T e is computed in the same way as in the model if fluxes are to be properly analysed.
For zonally-integrated vertical and meridional velocities W = xe xw w(x , θ, z)dx and V = xe xw v(x , θ, z)dx of a divergenceless flow we can introduce (see eg. Kundu et al., 2012) a streamfunction Ψ(z, θ) such that Here θ is the latitude in radians, z is the depth, R E is Earth's radius, v and w are the meridional and vertical velocities and 5 x w and x e the western and eastern boundaries in zonal direction. Following the definition, there are two convenient ways of computing global Ψ in geopotential coordinates: or

10
Here θ r the reference latitude. For global MOC computations it is the southernmost latitude of the Antarctic coast, where Ψ(z, θ r ) = 0. For regional MOCs, like the Atlantic MOC (AMOC), θ r is any convenient latitude where Ψ(z, θ r ) = 0 should be provided. For this, the equation 4 is usually used. In this equation the boundary condition is naturally taken into account by integrating from the bottom z = −H(x, θ). Note that both ways of computation are equivalent because the full velocity vector is divergence-free. In the following we discuss details of both methods of computation on unstructured meshes. Method

A (equation 3) involves vertical velocities and is more straightforward. Method B (equation 4) is based on horizontal velocities
and is slightly more complicated.

Method A
In FESOM2, the vertical velocity is conservatively remapped from vertices to cells using

20
where v(c) is the list of vertices of triangle c and N c is the number of the bottom level on triangle c. Indeed, it is easy to prove that v A kv w kv = c A kc w kc for FESOM2 discretization, so that the vertical (across level surface) transport is preserved.
Using triangles is more convenient in FESOM2 because bottom depth is constant on triangles. This remapping is not required in ICON-Ocean and MPAS-Ocean where the bottom depth is specified at scalar locations.
We introduce a set of latitude bins (θ i , θ i+1 ), θ i = θ 0 + i∆θ, i = 0, . . . , N θ covering the ocean domain. The computational 25 procedure is straightforward and is illustrated schematically in Fig. 2. triangle is in a bin if its centroid is in this bin. Triangles with centroids in dark blue, light blue and green fit in bins B1, B2 and B3 respectively.
-For each bin i find the list of triangles c(i) with centroids in these bins. They will be partly masked by bottom topography in deep layers, and we will formally write this list as c(ki), adding a layer index k. Subsequent computations are over triangles and levels, so that only c(ki) is needed.
-Compute ∆ψ ki as where c ki is the list of triangles the centers of which are in bin i at level k.
-Compute the meridional overturning streamfunction The procedure as written is strictly applicable in the case when level surfaces are fixed except for the surface. For z * vertical coordinates or for other options where where level surfaces are changing only slightly around their mean positions it can still 10 be used in most cases. It can be readily augmented with vertical remap to fixed levels by considering that the difference in transports (w kc − w (k+1)c )A c is linearly distributed within the layer in case when layers do not disappear, and level surfaces do not outcrop and stay at fixed depths where they cross topography. The method B should be used in more general cases.
Generally ∆θ should be taken about or larger than the typical size of triangles. The triangles that are counted as belonging to a bin are not necessarily confined to this bin, and the total area occupied by them differs from the bin area. However, 15 there generally are sufficiently many triangles in each bin, and one gets a smooth Ψ kj despite these effects. The procedure can be improved by conservative remapping into bins, which might be needed on coarse meshes. One may always check the bin attribution effect by repeating computations with smaller ∆θ. We also note that for instantaneous vertical velocities the procedure may result in Ψ different from zero at the surface. It will become zero only upon sufficient averaging which removes transient behavior of the surface.

5
The computations presented here can be generalized to some other sets of binning. Any sufficiently smooth scalar quantity defined at vertices or triangles can be used to introduce a set of bins. For example, being limited to the NA subpolar gyre, one may ask where the AMOC is forming using bins in mean sea surface height or barotropic streamfunction (see, e.g., Katsman et al. (2018)).
In the following we present an example showing differences between computations using different bins in ∆θ. For this, 10 FESOM was configured on a mesh with resolution varying from nominal one degree in the interior of the ocean to~1/3 degree in the equatorial belt and~24 km north of 50°N. We run the model for one year starting from climatology and compute the MOC from the annually averaged velocity. Because of starting the model at rest and short period of averaging we expect ∂η/∂t = 0, where η is the sea surface height. This, however, shall not affect the presented results. Using the coarsest bin size of 4 • the difference in MOC reaches locally above 5 Sv. As one would expect, decreasing the size of bins leads to convergence towards the solution obtained with the finest bin size of ∆θ = 0.125 • . We see that using bins of ∆θ = 0.25 • is already sufficient in this case because the mesh contains only few triangles that are smaller than the bin size.

Method B
Here the horizontal velocities are used. We select a set of latitudes θ i . The steps of the procedure are as follows.
-For each i draw a line θ = θ i and find a set of edges crossed by this line, as shown schematically in Fig. 6. For this, cycle through all edges, picking up those that satisfy the condition (θ v1 −θ i )(θ v2 −θ i ) < 0), with θ v1 and θ v2 the latitude of edge vertices. To avoid situations when the line passes exactly through the mesh vertex, a random noise of small 25 amplitude is added to the original θ i before edge e with vertices (v 1 , v 2 ) is tested. Schematic in Fig. 6 shows that the actual line through which transport is computed is a broken line composed of vectors (d ec1 , d ec2 ) related to the crossed edges. For a triangular C grid discretization one will deal with transports directly through the edges. The caveat in this case is that some of the crossed edges will be hanging and not contributing to the broken line. They are excluded by noticing that they have vertices that are encountered only once in the union of vertices of crossed edges. On hexagonal 30 C grids the procedure needs to deal with edges of dual triangular mesh. We denote the set of edges forming the broken line around θ = θ i as e(i). -The flux associated to the edge is given by the expression for F e above. The question now is the orientation of edges.
This question is trivially solved for each e by taking F e if θ v1 − θ i > 0 and −F e otherwise. It corresponds to keeping the normals to segments oriented so that transports are from the "northern" side of the broken curve. On triangular C grid the edge normal vectors used to introduce edge velocities can be selected as turned 90 • in positive direction from the edge direction. This will allow to solve the orientation problem in the same way.

5
-Since each of segments (d ec1 , d ec2 ) belongs to a particular cell, vertical integration is trivial for fixed level surfaces. If level surfaces are moving, the fluxes (transports) through the faces associated with segments are conservatively interpolated to the desired system of levels assuming linear distribution within model layers. In particular, the new system of levels can be specified in terms of potential density, with the result being the streamfunction in density coordinates. For each level the contributions from edges e ∈ e(i) are summed to get streamfunction at this level and the latitude θ i .

10
Note that the set of intersected edges may be ordered arbitrarily, the computation relies on the orientation of edges with respect to lines θ = θ i . This is the reason why the search for intersected edges remains relatively fast even on very large meshes.
Furthermore, it needs to be done only once for a particular mesh. Similarly to Method A, computations can be generalized to any set of lines, in particular to isolines of mean sea surface height or barotropic streamfunction. In both Methods we introduce masks if computations need to be confined to a particular basin.

15
Using this method we computed the streamfunction using the discrete spacing of ∆θ = 0.125 • . The difference to the streamfunction computed by method A is illustrated in Fig. 7. The discrepancy between both methods is caused by the difference of attribution of ocean volume to θ i . This, as shown in Fig. 7, can lead to a differences exceeding locally 1 Sv. These differences are not the errors, but uncertainty in the interpretation (see further).
If the modeled fluxes have been remapped onto the desired set of vertical levels as, for instance prescribed density lev-5 els, method B can be directly used for computing the MOC for a new vertical coordinate system. Figure  Note, that diagnostic of MOC in density coordinate can be also made in the same manner as method A. For this, the horizontal divergence needs to be remapped conservatively into density bins. From the horizontal divergence we then (1) diagnose the diapycnal velocity and (2) use it in method A.

Barotropic streamfunction
As follows from the equation for elevation, time mean vertically integrated horizontal velocity U is divergence free, ∇ η −H udz = ∇U = 0, i.e. it can be written in terms of the barotropic streamfunction Ψ as This streamfunction gives vertically integrated transport between two points at the surface.

Computations through binning
The barotropic streamfunction is more difficult to compute because binning has to be done in two directions. We introduce first a set of lines φ = φ j , where φ is the longitude, and φ j is the set of equally spaced longitude values over the basin of interest. As a first step the set of broken lines associated to each straight line φ = φ j is found. As the next step vertically integrated transports associated with the segments of broken line are computed. The final step is further binning of edges and 30 associated transports into equally spaced latitude intervals (θ i , θ i+1 ). Transport (and hence streamfunction) at each bin can be then computed by summing contributions going from the southern boundary where Ψ is set to zero.
This procedure can potentially be more noisy than computations of MOC, and may benefit from a conservative remap of the contributions from the segments in the second binning step (the number of segments in final bins is not necessarily large, in contrast to computations of meridional overturning).

5
According to the above procedure we computed the barotropic streamfunction using ∆θ, ∆φ = 0.25 • . Considering, that the procedure requires two-fold loop for (∆θ i , ∆φ j ) in case of large meshes and small bins it can become computationally heavy.
The result is illustrated in Fig. 8 (upper panel) and depicts reasonable structure of the main gyres with transports of 160 Sv and 70 Sv across Antarctic Circumpolar Current (ACC) and Gulf Stream, respectively.
In Fig. 8, middle and bottom panels show the differences between the streamfunctions if bins of 2 • and 1 • , respectively, are 10 used. As expected, the largest differences occur along the main gradients and reach of above 5 Sv along the ACC front. As in case with the MOC we note that these differences are not the errors, but uncertainty created by arbitrariness in the selection of bin size.

Computations through velocity curl
FESOM2 as its predecessor use implicit time stepping for the internal mode. The already available solver and routines need to 15 be only slightly adjusted to compute the barotropic streamfunction Ψ in the case when no-slip boundary conditions are applied.
In FESOM the discrete ζ is located at scalar points (at vertices), so modifications of the sea surface height solver to solve the above equations are indeed elementary. The difficulty in formal application of this approach is that the equation above needs 20 to be solved in a multiply connected domain with the Dirichlet boundary conditions provided on the periphery of each island and continent. Although these conditions can be formally provided by drawing lines connecting the islands and computing transports through the associated broken lines, this is tedious enough, especially when mesh resolution is high (and there are many islands). In the case of no-slip boundary conditions circulations along each island are identically zeros, and the equation above can be formally solved with the Dirichlet boundary condition on the southern boundary and the von Neumann boundary 25 condition ∂Ψ/∂n = 0 (n is the normal to the boundary). Although this condition does not ensure that Ψ = const over the periphery of any island, our experience with FESOM1.4 is that it works fine enough for practical purposes.
If we integrate the equation above over a scalar control volume (in FESOM2 scalar points are natural locations for relative vorticity ζ and streamfunction), we get 30 The contributions from edges on boundaries here are one-sided, including only segments that are wet (the first in the list in the case of FESOM). This automatically takes into account that there are no contributions from the boundary, as is the case for the no-slip boundary conditions. The operator on the left hand side in the case of FESOM is, up to the absence of depth weighting, the same as the part of operator used to compute the elevation, so the implementations is straightforward in the code (less so for post-processing). A clear drawback of this procedure is that it is not applicable for partial slip boundary conditions (it can be generalized, but will become too complicated). Since the methods based on bins was found to perform reliably, the curl-based method presents largely a historical interest. Computation of the barotropic streamfunction is currently implemented offline and from our experience it is slow because 25 of the loops along vertical and zonal directions are required. Hence we plan to implement the computation of the barotropic streamfunction following the philosophy of the in-situ computations (see e.g. Woodring et al., 2015).

Discussion
The general idea of the simple procedures described above is the use of transports as they are defined in an unstructured-mesh model, avoiding interpolation from an unstructured to a structured mesh. The diagnosed quantities such as meridional and 30 barotropic streamfunctions rely on the continuity equation, which is satisfied by the model only in a certain discrete sense.
Interpolation destroys this sense, requiring corrections and introducing interpretation errors related to these corrections. In practice the interpretation errors are significant, being on the level of Sverdrups for the meridional overturning as illustrated in Sidorenko et al. (2009), hampering discussions of MOC variability.
The algorithms above rely only on transports as they defined in models, and use conservative interpolation only in the vertical direction if required by a specified system of levels. We emphasize that the algorithms described are still sensitive to 5 parameter choices and thus contain interpretation uncertainty. In each case there is some sensitivity to how bins or vertical levels are selected. In Method B the straight line θ = θ i can be considered as centered in the respective bin, however the broken line drawn around the straight line is not necessarily centered within a bin. Drawing other possible broken lines in the bin is generally possible and can be proposed as a method to estimate this uncertainty. However, we would argue that such uncertainty is intrinsic to the diagnostics we are willing to compute. The true computation must rely on transport strictly consistent with 10 model discretization to avoid errors, and such transports are defined at irregular locations that generally do not lie on lines of latitude or longitude. A set of bins proposes some interpretation of integrated transports that is free of horizontal interpolation.
Any attempt to interpolate may create new uncertainties instead of making the analysis more accurate. These 'attribution' uncertainties have to be kept in mind especially in situations where small variability of MOC is the subject of analysis. Our experience thus far with the methods described above is that the computed patterns of MOC and barotropic streamfunction are 15 sufficiently smooth.

Conclusions
We describe a set of simple procedures intended to diagnose the meridional overturning and barotropic streamfunctions intended for unstructured meshes and requiring no interpolation of model output to regular meshes. We give application examples and discuss uncertainties involved. The procedures are described for FESOM2, but their adaptation for other discretizations 20 (MPAS or ICON) is straightforward. Our experience with using them indicates that they create much less difficulties with interpretation of model results than all our previous approaches based on interpolation.
Code availability. The code of the FESOM 2.0 model which was used to conduct the simulations for this paper is available at Zenodo (Sidorenko et al., 2020). The latest version of FESOM2 code can be downloaded from the public GitHub repository at https://github.com/FESOM/fesom2 under the GNU General Public License (GPLv2.0).

25
Data availability. Dataset related to this article can be found at Zenodo  Author contributions. DS and SD proposed the methods. DS and NK worked on the implementation, and NK, PS and QW contributed to the analyses. DS and SD wrote the the initial manuscript and all authors contributed to its final version.