Extracting Spatial Spectra Using Coarse‐Graining Based On Implicit Filters

Scale analysis based on coarse‐graining has been proposed recently as an alternative to Fourier analysis. It requires interpolation to a regular mesh for data from unstructured‐mesh models. We propose an alternative coarse‐graining method which relies on implicit filters using powers of discrete Laplacians. This method can work on arbitrary (structured or unstructured) meshes and is applicable to the direct output of unstructured‐mesh models. Illustrations and detailed discussions are provided for discrete fields placed at vertices of triangular meshes. The case with placement on triangles is also briefly discussed.


Introduction
The analysis of energy distribution over spatial scales, as well as the related analyses for energy generation, dissipation and transfer between scales are canonical tools in the study of eddy-driven flows.Typically, such analyses are conducted with the help of the Fourier transform.However, recently, an approach based on convolution, or coarsening has been proposed (Aluie et al., 2018;Sadek & Aluie, 2018).It not only provides energy distribution across scales but also gives the spatial distribution of energy at different scales.Additionally, this approach can be applied to flow domains of arbitrary shapes.Although the number of applications using this approach is still relatively limited compared to those utilizing Fourier spectra, they have demonstrated its usefulness.Some new-generation ocean circulation models such as FESOM2 (Danilov et al., 2017), MPAS-O (Ringler et al., 2013) and ICON-o (Korn, 2017) are formulated on unstructured triangular meshes, or their dual, quasihexagonal meshes.Scale or spectral analysis of fields simulated on such meshes typically requires interpolation to regular quadrilateral grids.Similar difficulties arise even with models using regular meshes on the sphere, as these meshes are rarely uniform, except in case where a relatively small region is explored and mesh nonuniformity can be ignored.local mesh inhomogeneity may pose challenges.Furthermore, the discrete Laplace operators can be implemented in full spherical geometry as used in models.For vector-valued data (horizontal velocities), vector Laplacians can be used either in a vector-invariant form or with metric terms, as will be discussed in some detail in this study.Our main intention below is to demonstrate the general utility of the method for unstructured meshes.
To start with, we explain the method using the one-dimensional case (Section 2).The extension necessary for horizontal fields is given in Section 3. The detail depends on the placement of discrete data and whether one deals with scalars or vectors.While our main focus is on the placement in Finite-volumE Sea ice-Ocean Model (FESOM), we also briefly discuss other placements.The analysis is applied to synthetic two-dimensional fields to illustrate technical aspects.The results are also compared to those obtained with the coarsening method proposed by Aluie et al. (2018).The final sections of the paper present discussions and conclusions.

Elementary Theory
The consideration presented here has very much in common with that in Guedot et al. (2015) and Sadek and Aluie (2018), because it incorporates elements from both approaches.Let ϕ(x) be a scalar field, with x lying in some domain D. We are interested in the distribution of the second moment (or variance if mean is removed) of this field over spatial scales.A coarse-graining approach, similar in spirit to that of Aluie et al. (2018) and Sadek and Aluie (2018) will be applied, but relying on implicit filters.The coarsened field ϕ ℓ (x) will be found by solving Here, Δ is the Laplacian, the smoothing scale is parameterized by ℓ, and γ is a parameter that tunes the relation of ℓ to wavenumbers, as will be explained further.The integer n defines the order of the implicit filter, which is second-order for n = 1 (in terminology of Guedot et al. (2015)), and order 2n in a general case.Equation 1 is complemented by boundary conditions on the boundary ∂D of the computational domain D, which will be discussed below.
A discrete Laplacian operator can be constructed independent of the type of the computational mesh and mesh geometry (flat or spherical), which makes the procedure suited for scale analysis on unstructured meshes.

1D Case
In order to clarify some details of the method, we begin with the one-dimensional case assuming periodic boundary conditions to facilitate the Fourier analysis.We assume n = 1 here, and explain later why larger n could be desirable.
In the Fourier space, where k is the wavenumber, the hat denotes Fourier transformed quantities and Ĝℓ (k) is the Fourier transform of the convolution kernel.The same notation k will be used for the absolute value of the wave vector in the discussion of energy spectra below in this section.Let E k be the Fourier energy spectrum of the field ϕ such that the mean energy over D is obtained by summation over positive k, where the integral is a replacement for the sum over Fourier modes.The area-mean energy of the coarse grained field, to be denoted by E ℓ , will be expressed as Clearly, E 0 = E.We can introduce (2) The quantity ɛ(ℓ) can be considered as energy distribution over scales ℓ because and ɛdℓ can be interpreted as the energy stored between the scales ℓ and ℓ + dℓ.This interpretation is not perfect because, in contrast to the Fourier spectrum, the contributions from different ℓ are not orthogonal, but it can still be useful in practice.
We will introduce one more derived quantity, proposed in Aluie et al. (2018) and Sadek and Aluie (2018).It can be more convenient as it presents a closer analogy to the Fourier spectrum.We define k ℓ = 1/ℓ and introduce The bar is used to distinguish this quantity from the Fourier energy spectrum.The form-factor F = ∂| Ĝℓ | 2 /∂k ℓ is a function of wavenumber k and has a peak at some k.It therefore picks up the energy E k in the vicinity of that peak wavenumber and ensures that E(k ℓ ) is an approximate representation of E k .It remains to be seen what is the correspondence between the peak wavenumber k and k ℓ , or the peak wavelength λ = 2π/k and ℓ, and when E(k ℓ ) represents E k sufficiently accurately.These questions will be addressed in the following sections.
As z = kℓ → 0, Ĝl (k) ≈ 1 z 2 / 24.We may require that all second-order kernels have such behavior in the vicinity of zero wavenumbers.For a Gaussian filter this requirement will lead to For the second order implicit filter this will lead to γ = 1/24, giving The left panel of Figure 1  √ ≈ 3.46 for the Gaussian and implicit second-order filter, and z ≈ 3.52 for the box filter.Taking the approximation k/k ℓ ≈ 3.5 we obtain ℓ = (3.5/2π)λ≈ 0.55λ, where λ = 2π/k is the wavelength.As expected, the scale ℓ is about half of a wavelength, and a more accurate statement has hardly a lot of sense since E k is a function of k, and the peak of the integrand in Equation 3 for given k ℓ will occur for somewhat different k than the peak of the form-factor.

Limitations
In order to see how E(k ℓ ) is related to the Fourier spectrum we assume that and zero otherwise.In this case, here z min = k min /k ℓ = k min ℓ and z max = k max /k ℓ = k max ℓ.The integrand in Equation 4does not depend on the scale ℓ, and if not for the variable integration limits, we would have concluded that E(k ℓ ) has the same spectral slope as E k just with a slightly different amplitude factor.However, the limits are functions of k ℓ .In practical applications with discrete data, k min = 2π/L x and k max = π/h, where L x and h are respectively the domain size and mesh cell size.
If γ 1/2 ℓ > h/π, which is generally the case, γ 1/2 z max > 1.Since the integrand in Equation 4 vanishes rapidly for γ 1/2 z > 1, the precise value of the upper limit will have only a small effect.The effect of the lower limit depends on the slope α via the numerator of the integrand.If α ≥ 2, the integrand takes the largest values at the lower limit, meaning that the integral in Equation 4 becomes sensitive to its lower limit, and that E(k ℓ ) will differ from E k for small k ℓ .The integral will be entirely determined by the behavior at the lower limit for α > 3.Such slopes cannot be diagnosed with second-order filters, which will return E(k ℓ ) ∼ k 3 ℓ .Sadek and Aluie (2018) give an illustration of such a saturation effect for slopes steeper than 3.
Note that even if the Fourier energy spectra are sufficiently flat (e.g., with α = 5/3), the integral in Equation 4 can still be sensitive to the lower limit because z 2 α varies too slowly.In a more general situation, when E k has a local spectral peak at intermediate scales, the position of a local peak in E(k ℓ ) could be aliased because the form-factor is rather broad for the second-order filters (see Figure 1, right panel).This is the reason why it is advisable to use higher orders n in Equation 1. Indeed, Equation 4 for arbitrary n transforms to

E(k
(1 + γz 2n ) 3 . (5) For n = 2, that is, for the biharmonic filter, one will need to take γ = 1/288 to have the same position of the peak in the form-factor as in the second-order case with γ = 1/24 because the peak is always at 2γz 2n = 1.This case is depicted in Figure 1 with the dashed line, and the right panel illustrates that the peak becomes sharper than for the second-order filters.Importantly, the form-factor remains closer to zero at small wavenumers and is well-suited for a practically interesting range of slopes α ≤ 3. The relation between the Fourier spectrum and the scale spectrum becomes more local.Sadek and Aluie (2018) discuss the construction of high-order explicit convolution kernels and illustrate the improvements in the performance of the scale analyses with high-order kernels.We present similar illustrations for the biharmonic case later.
Note that one is free in the selection of parameter γ.For practical purposes it might be convenient to take γ = 1/2, in which case the peak is at k = k ℓ and ℓ acquires the sense of inverse wavenumber, that is, ℓ = λ/(2π), where λ is the wavelength.
Although the form-factor is narrower for the biharmonic filter than for the harmonic one, it has a finite width as a function of k.The form-factor for the n-harmonic implicit filter is given by (1 + z 2n /2) Its width at the half amplitude level can be found from which gives z 2n ≈ 0.2 and z 2n = 4 for the left and right boundaries of the peak.The order of the filter determines the actual width.The peak is approximately between 0.45 ≤ z ≤ 2 for the Laplacian filter, 0.67 ≤ z ≤ 1.4 for the biharmonic filter, and 0.82 ≤ z ≤ 1.2 for the bi-biharmonic filter.For the harmonic filter the peak is wider than two octaves, and for the biharmonic filter it is still wider than one octave.The coarse-graining method would not be able to isolate peaks in the Fourier spectra if they are sufficiently narrow unless a sufficiently high-order filter is used.An illustration is given below.For completeness we note that the width of the main peak of box filter formfactor is about 0.7 of the width for Laplacian filter.It is still 1.5 times wider that the width for the biharmonic filter.

Modifications in the Discrete Case
The eigenvalue of the discrete Laplacian is an approximation to k 2 of the continuous case.In the onedimensional case the simplest numerical implementation of the Laplacian operator on a regular mesh (Lϕ) j = ϕ j 1 + ϕ j+1 2ϕ j )/ h 2 , where j is the cell index and h the cell size, leads to the Fourier symbol L k = (4/ h 2 ) sin 2 (kh/2).For the largest resolvable wavenumber k = π/h, the ratio of discrete to continuous Fourier symbols is L k /k 2 = (4/h 2 ) sin 2 (π/2) (h/π) 2 = 4/π 2 , which is not negligible.
We explore the consequences of this difference, by writing where the scale parameter β = β(k ℓ ) is introduced formally instead of ℓ 2 = 1/ k 2 ℓ .In this case the form-factor is Considering the form-factor as a function of wavenumber k, we find that it reaches maximum at Journal of Advances in Modeling Earth Systems We take γ = 1/2 everywhere further to simplify this condition.
There are two obvious options.We can modify the scale parameter β taking β = β M = 1/ L k ℓ , which is similar to the correction used by Guedot et al. (2015).In this case the form-factor peaks at k = k ℓ = 1/ℓ, so that k ℓ acquires the sense of the standard wavenumber, which lies in the range from 2π/L x to π/h.The drawback of this choice is that β 1 ∂β/∂k ℓ ∼ cos(k ℓ h/2)/sin(k ℓ h/2), which tends to zero at k ℓ = π/h, artificially damping the scale spectrum near the cutoff wavenumber π/h.
The other choice is to leave β = 1/ k 2 ℓ .In this case the position of the peak of the form-factor is at L k = k 2 ℓ .The relationship between k and k ℓ depends on the eigenvalue of discrete Laplacian.For the standard 1D Laplacian we find sin(kh/2) = k ℓ h/2.The peak at k = π/h corresponds to k ℓ = 2/h and the peak at k = π/(2h) corresponds to . The latter values are already close, and the correspondence between the peak k and k ℓ is recovered for smaller wavenumbers.
As is seen, both choices introduce certain complications at grid scales.However, discrete operators in ocean circulation models also contain errors at grid scales, so that the scale spectra at such scales are of limited interest.If one considers wavenumbers smaller than π/(2h), both choices of β and respective scale spectra become close.One can also ignore the distinction between k ℓ and the peak k in this range if β = 1/ k 2 ℓ , and this choice of β becomes a simpler option.We will, nevertheless, keep the option β = β M below for completeness.Guedot et al. (2015) use the 1D Laplacian eigenvalue above also on unstructured meshes, identifying h with the length of the side of a triangle, but we opt for a slightly different value, see below.Since h is variable on unstructured meshes, β M (if used) should be included between divergence and gradient in the scalar Laplacian, and similarly in the vector Laplacian.
For the comparison with the box filter, we recall that the peak of the form-factor is at kℓ box ≈ 3.5, which means that ℓ box ≈ 3.5/k ℓ in terms of k ℓ used further for implicit filters.

Unstructured Meshes
Laplace operators on unstructured triangular meshes depend on the placement of discrete degrees of freedom.Some of them are reviewed by Danilov and Kutsenko (2019) and Klemp (2017).As applied to the existing unstructuredmesh global ocean circulation models, one needs to distinguish between the cases of vertex, cell and mid-edge placement, and the cases of scalar and vector operators.The main focus below will be on the FESOM discretization which places scalars at vertices and horizontal velocities on triangles.Although horizontal velocities in FESOM are placed on triangles, in the default implementation we seek coarse grained velocities at vertex locations because this leads to matrices with a twice smaller dimension.One may, nevertheless, need the native placement to analyze dissipation spectra, but we do not consider it in detail here, limiting ourselves only to the case of scalars placed on triangles.The particular feature of this case is that the Laplacians either contain spurious eigenvalues or a non-trivial kernel.We mention that the cell placement of scalars on hexagonal meshes in MPAS-Ocean differs from the vertex placement of FESOM only in the detail of control volumes, and that scalars are placed on triangles in ICON-o.Normal component of velocities are placed at mid-edges in MPAS-Ocean and ICON-o.In these cases vector-invariant forms of vector Laplacians are appropriate, which will not be analyzed in detail here.
The discretization of Laplacians on unstructured meshes can be obtained by finite volume (as in Guedot et al. (2015)) and finite element methods, and we briefly describe both for the vertex case.The discretization of vector Laplacian for FESOM which works with the full horizontal velocity vector (quasi B grid) can be done similarly to scalars, but with the addition of metric terms as discussed below.Using the vector invariant form of the vector Laplacian will not lead to advantages for the quasi B grid.

Scalars, Vertex Placement
One can use either the finite volume or finite element method for discretization of the Laplacian operator.They lead to slightly different results on non-uniform meshes.In the first case, which is also the choice in Guedot et al. (2015), the minus Laplacian operator D, modified with the scale factor β, is where the index v enumerates vertices, N(v) is the set of neighboring vertices, l vv′ is the length of the edge between vertices, n vv′ is the unit vector from v to v′, s vv′ is the boundary segment vector and A v the area of scalar mediandual control volume containing vertex v, see Figure 2. We included the scale factor β vv′ as an edge-based quantity.
It can be either β M vv′ (see below) or 1/ k 2 ℓ , in which case it can be taken as a factor before the sum.The quantity s vv′ is a combination of two contributions obtained by rotating vectors connecting mid-edge to cell centroids in the case of FESOM (median-dual control volume), as illustrated in Figure 2, but it would be related to a single segment on Voronoi (hexagonal) meshes.Because of A v on the left hand side, D is not symmetric (with respect to v and v′), but the right hand side of Equation 7gives a symmetric operator.This is why Equation 1 has to be discretized in area-weighted sense (see below).
Assuming ϕ ℓ ,ϕ ∼ e ik ⋅ x , where k = (k, l) is the wave vector and x the coordinate vector, and further assuming that the mesh consists of equilateral triangles with the side a and height h = a ̅̅ ̅ 3

√
/ 2, with one side oriented along the xdirection, we obtain the Fourier symbol of the vertex Laplacian ( D for β vv′ = 1) The Fourier symbol depends on the wave vector orientation.The boundary of the first Brillouin zone is closest to k = 0 in six directions (one is k = 0), and in such directions L k = L K = (4/h 2 ) sin 2 (Kh/2), where K = |k| (see Danilov, 2022).This expression is similar to the one discussed for the one-dimensional case, with the difference that h is not the length of the triangle side, but the triangle height.We take β M vv′ = L k ℓ , where L k ℓ = L K | K=k ℓ , and express the height assuming regular meshes, h = ̅̅ ̅ 3 √ l vv′ / 2. The wavenumber k ℓ is the absolute value of twodimensional wave vector and should be taken between 2π/L x and π/h.We will suppress the subscript ℓ in ϕ ℓ further for compactness.The discrete version of Equation 1 for n = 1 becomes where summation over the repeating index v′ is implied.The system matrix is Figure 2. Finite-volume computations of Laplacians for the vertex (left) and cell (right) quantities.In the vertex case, n vv′ is the unit vector in the direction of edge, and s vv′ = s l + s r , where s l and s r are the left and right vectors from the mid-edge to cell centers rotated to outer normal directions.In the cell case, n cc′ is the unit vector along the line connecting cell centers and s cc′ is the edge vector rotated toward the outer normal.
with I vv′ the identity matrix and D vv′ the matrix of the operator in Equation 7.
Taken in this form, S vv′ is a positive symmetric matrix, and solution for ϕ can be obtained by the conjugate gradient method.
At this place we turn to the question of boundary conditions.The procedure of constructing the Laplacian, described above, implies that the component of gradients that are normal to domain boundaries are thought to be zero, that is, we are dealing with the zero Neumann boundary conditions.One can consider the Dirichlet boundary conditions, but requiring that ϕ = ϕ at boundary points will introduce small scales at the boundary and break the symmetry of the system matrix.One more choice is to extend the mesh into land and continue ϕ there (i.e., through the nearest interpolation), which might reduce the sensitivity to boundary conditions.These possibilities are not explored here.We use the zero normal flux boundary conditions in all cases below to ensure the symmetry and positivity of system matrices.Such systems of linear equations can be solved with the conjugate gradient method.
In the case of a biharmonic filter the matrix of the biharmonic operator can be obtained by applying the procedure used for D twice, with summation over v′ implied.Although D v′v ′′ includes division over A v′ , the combination A v D vv′ D v′v ′′ is symmetric.This procedure also implies the additional boundary condition of zero higher-order normal flux n ⋅ ∇Δϕ ℓ = 0 on the domain boundary.Here n is the unit outer normal.The higher-order operators n > 2 can be obtained similarly, with the additional boundary conditions of zero normal high-order fluxes.While the area weighting in the matrices above is important for symmetry, it will be eliminated by diagonal preconditioners, meaning that S vv′ can be assembled without A v in this case.
Equation 1 can also be discretized with the help of the finite element method.Although this method leads to a similar discretization for scalars as the finite volume method, it is our preferred approach because it is more convenient for vector Laplacians as concerns the metric terms.Consider first the case n = 1.We obtain the weak formulation of Equation 1 by multiplying it with some sufficiently smooth function w(x), integrating over the domain D and transforming the Laplacian term by integration by parts, here the boundary term appearing after the integration by parts is set to zero, which corresponds to natural (Neumann) boundary conditions, similarly to the finite-volume derivation.We write the discrete fields as series ϕ = ∑ v′ ϕ v′ N v′ (x) and ϕ = ∑ v′ ϕ v′ N v′ (x), where N v (x) is a piecewise linear basis function.This function equals 1 at the position of vertex v, decays to 0 at vertices connected to v by edges, and is zero outside the stencil of nearest triangles.To obtain the continuous Galerkin approximation we require that the weak equation above is valid if w is any of the N v , where summation over repeating v′ is implied, M vv′ = ∫ N v N v′ dS is the mass matrix and Dvv′ = ∫ β∇N v ⋅ ∇N v′ dS.

The diagonally lumped approximation of mass matrix is M
corresponds to D and we get the same final discretization on uniform equilateral meshes as for the finite volume method except for the treatment of β.Since computation of integrals are done on triangles (mesh cells), the scaling parameter is treated as a cell quantity if the modified option is used.Gradients of linear basis functions are constant on triangles, and their integration is just multiplication with triangle area.The height of triangle in β M on cell c is conveniently expressed in terms of triangle area as h = 3 1/ 4 ̅̅̅̅̅ A c √ .Once again, we have symmetric matrices.

The implementation of biharmonic operator leads to
Here, in reality, the procedure consists of two steps.We explain them ignoring β for brevity.We express ΔΔϕ as Δψ, where ψ = Δϕ.Discretizing the latter with the finite element method we get M L vv′ ψ v′ = Dvv′ ϕ v′ .In the second step we discretize ϕ + (1/ 2)Δψ = ϕ, which leads to the biharmonic filter written above.On both steps the flux terms are omitted, which is equivalent to the boundary condition of no normal flux (first step) and no normal high-order flux (second step) at the domain boundary.The same as with the finite-volume method, the procedure can be generalized to higher-order filters.

Spherical Geometry and Vertex-Based Vector Laplacians
FESOM uses longitude-latitude coordinates in spherical geometry (for detail, see Danilov et al., 2017).The cosine of latitude is approximated by a constant on triangles and local distances on each triangle are calculated in a local Cartesian frame, with x and y aligned with local zonal and meridional directions.The areas of triangles are computed assuming that triangles are flat, and areas of scalar control volumes are obtained by taking 1/3 of the area of all triangles joining at a given vertex.This way of treating spherical geometry is fully compatible with the finite element computations described in the previous section.For the finite volume treatment, the length of edges is determined using the cosine at mid-edges.
Although horizontal (zonal and meridional) velocities in FESOM are placed at triangles, in the present implementation the filtered velocities are calculated at vertices because there are twice fewer vertices than triangles, and matrices to be iteratively inverted have twice lower dimension.
Let u = (u,v) be the filtered horizontal velocity, and u = (u, v) the original velocity (we suppress the ℓ subscript here too).We will use index v for vertex-based velocities, and index c for cell-based original velocities.The weak equation will take the form where now w is the vector-valued smooth function.Writing u = ∑ v′ u v′ N v′ (x), u = ∑ c u c M c (x), where M c is the indicator function equal 1 on triangle c and 0 outside it, and requiring that the weak equation is valid for w = w v N v for any vertex v, we obtain the set of discrete equations.However, because the unit zonal and meridional vectors are varying, there will be additional metric terms in the expression for velocity gradient, and similarly in the expression for the gradient of test function w.Here, m = tanθ/R e , with θ the latitude and R e earth's radius.Combining all metric terms, we will get the following x and y equations where Dvv′ now includes an additional term due to metric, ) dS is the operator accounting for metric terms linear in m and R vc = ∫ N v M c dS is the right hand side operator which is the projection operator onto the space of piecewise linear functions.Summation is implied over repeating indices in matrix-vector products.
The metric terms couple equations for u and v, and the dimension of system matrix becomes twice larger than for scalar fields.While it can still be solved, the metric terms can be accounted approximately if wavelengths 2π/k ℓ are not very large by putting them on the right hand side and estimating them from the solution for the nearest available k ℓ .For the largest k ℓ they can be estimated using the original field.On basin scales, this approximation is not necessarily sufficient and the full system has to be solved.For higher-order filters, we introduce ).
The system matrix in the biharmonic case becomes

Scalars, Cell Placement
The smallest-stencil negative Laplacian D in this case is obtained using the finite volume method and is written similarly to Equation 7 Here summation is over the neighbors of cell (triangle) c, which are the cells that have common edges with c (see Figure 2), l cc′ is the distance between the cell centers, A c is the cell area, n cc′ is the unit vector in the direction from c to c′ and s cc′ is the cell side vector rotated to the direction of the outer normal.For meshes used in discretizations of C-grid type, the c points are associated with circumcenters.In this case the vectors n cc′ and s cc′ have the same direction, and this expression presents the divergence of gradient taken with minus sign.If c is associated with cell centroids, on distorted meshes one accounts only for the component of gradient in the direction of n cc′ while full gradient is formally needed (discretization given by Equation 7 has similar limitations).We have found that in this case a simplified version of Equation 10 where A = ̅̅ ̅ 3 √ for equilateral triangles and 3/2 for rectangular ones leads to a better convergence.More precise computations are of course possible, but they will not be considered in detail here.Although details of discretization introduce some uncertainty in the analysis, they will basically affect only grid scales.
To estimate the eigenvalues of D, we use a uniform equilateral mesh, setting β cc′ = 1 in Equation 11.In order to find the eigenvalues we first note that for any two neighboring triangles the orientation of stencils of neighbors depends on the orientation of triangles themselves.For this reason, an elementary Fourier harmonic is described by two amplitudes and the Fourier symbol becomes a 2 by 2 matrix (see e.g., Danilov & Kutsenko, 2019;Klemp, 2017), with two eigenvalues The plus sign corresponds to the physical eigenvalue which tends to k 2 l 2 in the limit of small wavenumbers.The minus sign corresponds to a numerical eigenvalue which is anomalously high in absolute value (it tends to 24/a 2 in the limit of small wavenumbers).
In the direction k = 0, |V| 2 = 5 + 4 cos(Kh), which corresponds to a more accurate behavior than in the vertex case above, giving the physical L k equal to 6/h 2 at the largest wavenumber which is closer to π 2 /h 2 of the continuous case.The eigenvector of the numerical eigenvalue corresponds to a grid-scale pattern that may potentially be present because of a slight difference in differential operators on differently oriented triangles.However, even if data contain some contributions that project on the numerical mode, these contributions will be strongly damped in the filtered fields if the implicit filter is used.In contrast, the presence of spurious eigenvalues could be problematic for explicit filters based on series in Laplacians.
An alternative approach is to use a discrete Laplacian that is based on a wider stencil.An obvious choice is to reconstruct full gradient vectors on triangles based on the values on three neighboring triangles.Gradients are then averaged to edges, and their divergence gives a Laplacian (see, e.g., Danilov & Kutsenko, 2019).It turns out that such a discrete Laplacian has the same physical eigenvalue as the vertex Laplacian above, that is, it is less accurate.Its numerical eigenvalue is zero.This Laplacian would be less useful for implicit filter applications because it will preserve grid-scale noise in coarse-grained fields.A combination of this Laplacian with the smallest-stencil Laplacian can be a way to reduce spurious eigenvalue for applications with explicit filters.

Other Discretizations
The scalar operator in Equation 10 can be applied to the components of cell horizontal velocity in plane geometry.
Metric terms need to be added in spherical geometry, we do not discuss them here.The vector-invariant form of the vector Laplacian, Δu = ∇∇ ⋅ u curl curlu would allow one to avoid the question on metric terms.It depends on computations of the divergence and inner curl.For cell velocities a natural way would be to compute the divergence and curl at vertex control volumes using the divergence and circulation theorems because velocities are known at boundaries of such control volumes.This would lead to an analog of the wider-stencil Laplacian discussed in the previous section, which is not necessarily optimal because of its large kernel.
C-grid types of discretizations used in MPAS-O and ICON-o will rely on the vector-invariant forms of vector Laplacians.These Laplacians have large spurious eigenvalues (see e.g., Danilov & Kutsenko, 2019;Klemp, 2017), but for the reasons mentioned above, this is rather an advantage.

Comments on Solving Procedure
As recommended by Guedot et al. (2015), instead of solving where ϕ and ϕ are the vectors of discrete values, one works with perturbations ϕ′ = ϕ ϕ, with the result Sϕ′ = ϕ Sϕ.
In this case the right hand side has emphasis on small scales, and ϕ′ is the fine-scale contribution that has to be removed from the original data.This is expected to improve the convergence, which is more difficult for large scales.It indeed leads to some improvements in our case, however the main issue is that conditioning of S worsens as k ℓ decreases.
Also, as explained by Guedot et al. (2015), one will split the operators when n > 2. For example, and and so on (γ = 1/2).In principle, using tri-harmonic and bi-biharmonic operators only doubles the work of iterative solvers once matrices of harmonic and biharmonic operators D and D 2 are available.While biharmonic operator should be sufficient in most cases for general conclusion about energy spectra, one may opt for higherorder in the case when precise position of spectral peaks has to be estimated.

Illustrations
For illustrations here we first generate a scalar field ϕ(x) on a regular quadrilateral mesh covering a square domain with size L x = 1,024 km.The initial mesh has the resolution a = 4 km (the side of quadrilateral cell).The field ϕ is initially specified at vertices of this mesh by taking random values uniformly distributed around zero.It is Fourier transformed, and each wave component is divided by K (α+1)/2 , where K is the absolute value of the wave vector.After that, the inverse Fourier transform is carried out.In this way, fields with Fourier power spectra E K ∼ K α are created, and α = 2 is used below.The Fourier amplitudes are additionally modified to add peaks in Fourier power spectra in Section 4.2.
One triangular mesh is obtained by bisecting the quadrilateral cells of the primary mesh.The other one is an equilateral mesh with triangle side a.The third and the fourth meshes are unstructured with the triangle sides varying between a and 3a and a and 5a respectively.These meshes will be denoted further as Q, E, U1 and U2.The length of triangle side of U1 and U2 is shown in Figure 3.The smallest triangles (with the side of 4 km) cover about 25% of the area on these meshes.Two circular regions are included to modify the boundary of U1 and U2.One more unstructured mesh with the triangle sides varying between a and 2a, but covering a circular area, is used for tests in spherical geometry.The largest resolvable wavenumber is defined by the small side of triangles on mesh Q and the height of triangle on mesh E.Although the cells of U1 and U2 vary, spectra are computed up to the same wavenumber as on mesh E. The initial field on mesh Q is the constructed one, and it is linearly interpolated to meshes E, U1 and U2, respectively.Except for Figures 4 and 9, we use the modified β.We take γ = 1/2 which ensures that k ℓ corresponds to the ordinary wavenumber and ℓ has the sense of inverse wavenumber.A conjugate gradient solver is applied in all cases.In all figures showing spectral density, the vertical axis is in units of the field squared times km/cycle, and wavenumbers on the horizontal axis are in cycle/km.

Uniform and Distorted Triangular Meshes
Two panels of Figure 4 present the computed scale spectra on triangular meshes Q, E, U1 and U2 using harmonic (left) and bi-harmonic (right) filters for realizations of field ϕ with the power density slope α = 2. Here, the scaling factor β = 1/ k 2 ℓ .The Fourier spectrum is in red, it contains fluctuations because it is the spectrum of a realization.The dashed black line shows the slope 2. The same realization of ϕ is used in the left and right panels.The blue, solid black, gray and green lines give scale spectra on meshes Q, E, U1 and U2 respectively.The horizontal axis correspond to K for the Fourier spectrum and k ℓ for the scale spectra.The k ℓ values are sampled up to π/h, where h is the smaller side of triangle on mesh Q or the height of the smallest triangle on other meshes.Only the wavenumbers smaller than the cutoff of the Fourier spectrum should be compared.Since β = 1/ k 2 ℓ , we should stop even earlier at k ℓ = 2/h, interpreting it as K = π/h.As expected, for the original data (blue curves) the bi-harmonic filter ensures a closer agreement between the scale and Fourier spectra.As already mentioned, the biharmonic filter can also be used for sufficiently steep spectra (formally α ≤ 5), while the use of harmonic filter is limited to α ≤ 3. We checked this, but do not demonstrate it here.
The scale spectra on meshes E, U1 and U2 (solid black, gray and green curves) deviate from the scale spectrum on mesh Q for large k ℓ .For mesh E this is the consequence of using interpolated data, which have a smaller variance on grid scales.The spectra for U1 and U2 meshes deviate stronger than for mesh E because fine triangles occupy only about a quarter of these meshes.The spectra start to fall off when k ℓ exceeds the cutoff defined by the  .Shown are the Fourier spectrum (red), the scale spectrum on a uniform triangular mesh obtained by bisecting quadrilateral cells (blue, mesh Q), on a regular equilateral mesh (solid black, mesh E) and on irregular unstructured meshes U1 and U2 (gray and green, respectively).The wavenumber is in cycle/km.The horizontal axis corresponds to K for the Fourier spectrum and to k ℓ = 1/ℓ for scale spectra.coarsest mesh triangles.The differences between the spectra are smaller for the harmonic filter.This is by all probability linked to its broader form-factor, which picks up larger contributions from smaller wavenumbers where the spectral density is higher.Since small triangles occupy only one quarter of meshes U1 and U2, this might lead to a four times smaller spectral energy density at the smallest scales.The observed difference with respect to the black curve of mesh E is smaller even for the biharmonic filter, which once again is due to the width of the form-factor.The difficulties with spectra on unstructured meshes can be avoided either by stopping at wavelengths resolved by the entire mesh, or by limiting the analysis to the high-resolution patch.
The scale spectra tend to bend down compared to the Fourier spectrum on the side of smallest wavenumbers, which is presumably because a part of formfactor curve is cut on the left.Similar tendency on the side of large wavenumbers is partly related to the difference between K and k ℓ in this case.
Apart from effects caused by interpolation or created by mesh nonuniformity, the scale spectra are insensitive to the mesh type, which proves that the method can work on different meshes, including unstructured meshes.The scale spectra are always smoother than the realization of the Fourier spectrum, as they correspond to the integration over some range of Fourier wavenumbers in the vicinity of k ℓ , which smooths the fluctuations of the Fourier spectrum.On a practical side this means that one may need less averaging over field realizations as compared to the Fourier spectrum.

Peak Detection
All further examples except for Section 4.5 use the modified β, so that the correspondence between the wavenumbers K and k ℓ is maintained.Fields with peaks in the Fourier spectrum have been constructed following the same procedure as above, but applying additional factors to the amplitudes of Fourier harmonics.Figure 5 displays two examples illustrating possible behavior.We use mesh Q and the biharmonic filter.In the left panel, the peak in the Fourier spectrum is accurately recovered in the scale spectrum.In contrast, the close peaks in the scale spectrum in the right panel are not detected, and the peak around k ℓ ≈ 0.2 cycle/km is nearly missed.The success in the first case is because the Fourier spectrum decays sufficiently fast on both sides of the peak and there are no other close peaks.The form-factor responsible for the correspondence between the Fourier and scale spectra is still not sharp enough for the field used in the right panel, even for the biharmonic filter, and higher-order filters have to be used if finding exact position of peaks is of utmost importance.The estimates of the width of the form-factor from Section 2.3 can be used to see that the high-wavenumber peaks in the Fourier spectrum are too narrow to be resolved.

Harmonic Filter Versus Explicit Box-Type Filter
We use a box-type filter with the kernel where A is the normalization constant selected so as to ensure that ∫G(x) dx = 1, a is the (fixed) mesh cell size, c = 1.3 is a numerical factor controlling the sharpness of the box filter transition, and 1.75/k ℓ is equivalent to the scale ℓ box /2 of the box filter, as discussed above.Figure 6 compares the scale spectrum obtained by explicit coarse-graining with this box-type filter (blue line), with the scale spectrum obtained by Laplacian smoothing (solid black line) on mesh Q.Both methods lead to nearly the same spectra, with the largest differences on the side of large wavenumbers.The selection of parameter c was experimental.Decreasing it increases the discrepancy   12) (blue) compared to scale spectrum computed with Laplacian smoothing (solid black).The Fourier spectrum is in red, and the dashed line corresponds to the slope of 2.
between the two curves, which is still smaller than the deviation from the Fourier spectrum.If c is increased, the box filter transition becomes sharper, but this leads to some unevenness of the spectrum on the side of largest wavenumbers.We conclude that both methods agree, as expected, and that the attribution ℓ box = 3.5/k ℓ works well.The excessive decay at the spectral end (k ℓ ≥ 0.5 cycle/km) shown by solid black line is the consequence of using the modified β in the analysis.

Spherical Geometry
An unstructured variable-resolution mesh covering a circular area with diameter 1,024 km is created.The sides of triangles vary between a = 4 km and 2a.The realization of ϕ(x) with 2 power spectrum is linearly interpolated on this mesh.The mesh coordinates are then approximately transformed to longitude and latitude by multiplying the coordinates with 180/(πR e ).The sphere is then rotated in the meridional direction by π/6, π/3, and 5π/12 to produce three variants of the same mesh having various longitude and latitude coordinates, as shown in Figure 7.We apply the coarse-graining based on the biharmonic filter to these three meshes.Although they correspond to the same original mesh and data, they have different coordinates on the sphere and the cosine of latitude varies much stronger on the higher-latitude mesh than on the other two meshes.However, the computed scale spectra, shown in Figure 8, are indistinguishable.This indicates that the approximation of cosine by a constant value on triangles is sufficient.By using a circular mesh in this example we illustrate that the smoothing procedure can work in a domain of arbitrary shape.

Kinetic Energy Spectra in Spherical Geometry
For this test case, a realization of scalar field ϕ(x) with power spectrum with a 4 slope was created on mesh Q.The mesh coordinates were transformed to longitude and latitude by multiplying them by 180/(πR e ).To make the metric terms larger, the mesh was moved along the zeroth meridian so that its south-west corner is at (0°, 75°).After this step the mesh becomes stretched zonally if viewed in longitude-latitude coordinates on the sphere.The scalar field ϕ(x) is interpreted as a streamfunction given at mesh vertices.Velocities are computed at the centers of triangles by calculating derivatives of the streamfunction given at vertices.We solved the full system of Equations 8 and 9, building matrices with twice the dimension of the scalar problem.Blue curves in Figure 9 show the scale spectrum of kinetic energy obtained with the harmonic filter (left) and biharmonic filter (right).The spectra computed by solving two problems for velocity components with metric terms ignored are hidden behind the blue curves and only begin to show up (in black) at the smallest wavenumbers.In the example used, solving Equations 8 and 9 takes nearly the same computational time as solving two smaller-size problems for velocity components separately, but requires more storage.Note that the computed kinetic energy spectrum should be steeper than 2 at large wavenumbers because numerical differentiation of the streamfunction adds K 2 to the streamfunction power spectrum only on relatively large scales.Additional attenuation is introduced by averaging on the rhs of Equations 8 and 9.These effects are only seen with the biharmonic filter.

Discussions
In this work we present an elementary theory of scale analysis based on implicit filters that use powers of Laplacian.Since a discrete Laplacian can be written for any mesh type, the procedure is sufficiently general.Our focus was on triangular meshes and on the scalar and vector data placement of FESOM,  but the definition such as Equation 7 are valid for arbitrary meshes.The use of longitude-latitude coordinates was dictated by their convenience for the discretization of FESOM.The method can be cast in the form that relies on local tangent-plane calculations, which is a subject of future work.Such a form will be more convenient for the discretizations of C-grid type (MPAS-Ocean and ICON-o) and can be applied even to fields on the entire sphere.
The formulation depends on the discretization of Laplacians and on the imposed boundary conditions.We assume that the normal fluxes are zero at the flow boundaries for the harmonic filter, and additionally require that the normal higher-order fluxes are zero for higher-order filters.This leads to symmetric and positive-definite system matrices.In this case the coarsegraining problem is solvable with the conjugate gradient method.More general boundary conditions will destroy the symmetry of system matrices and were not tried here for this reason.The question on whether they have some advantages is a topic for further studies.While no-flux boundary conditions agree with the boundary conditions imposed on tracers in ocean circulation models, this agreement plays no special role here as the equations solved in the models and in the implicit filters are different.A valuable side effect of these conditions is the conservation of total tracer content (total momentum or angular momentum in flat or spherical geometry, respectively) in coarsegrained fields (see also Grooms et al., 2021).
Note that the explicit coarse-graining also depends on the definition of convolution in the vicinity of domain boundaries, so the dependence on the boundary conditions is common aspect of coarse-graining methods.
Discrete Laplacians deviate from their continuous counterparts on grid scales and may have spurious modes in a general case.We used only the simplest options for discrete Laplacians and explored in some detail the case of vertex placement.The deviation of eigenvalues of discrete Laplacians from the continuous case creates some difficulties in the interpretation of scale spectra near the spectral end.Our estimates of eigenvalues were based on the Fourier analysis, which can only give an approximate result because the eigenvalues of discrete Laplacians depend on the boundary conditions and mesh structure in each practical case.Our analysis of possible effects is therefore only approximate, but it gives an idea of how the spectra are affected.It suggests that using the modified scale factor similar to Guedot et al. (2015) one can keep the correspondence between k ℓ = 1/ℓ and the standard wavenumber K, but creates an artificial spectral decay, as seen in Figures 6 and 8.If the original definition β = 1/ k 2 ℓ = ℓ 2 is retained, the direct correspondence between K and k ℓ is violated near the spectral end.While these difficulties are not important in practical cases dealing with the behavior of well-resolved scales, one may be interested in studying energy pileup at grid scales, where the interpretation of spectral density might be ambiguous.Note that similar issues may also accompany the coarse-graining based on explicit box-type filters when ℓ box /2 becomes close to the width of transition part (see Equation 12).
While more accurate discrete Laplacians will reduce these difficulties, the construction of such Laplacians for general unstructured meshes is far from straightforward.However, one does not always need the exact correspondence between the Fourier and scale spectra: it may be sufficient just to compare the spectra in terms of k ℓ .
The interpretation of scale spectra on unstructured meshes requires some care because of mesh variability.If the mesh resolution varies substantially within the domain being studied, going to the scales smaller than the wavelength determined by the coarse cells of the mesh is not generally reliable, except when the coarse part occupies only a fraction of domain.Laplacian filters may mask the problem, as is illustrated in Figures 4 and 9, but this does not make the analysis more reliable.A better solution is to use a subdomain where the mesh is approximately uniform.Since the coarse-graining procedure produces ϕ ℓ at every location, the calculation of E(k ℓ ) can be easily constrained to a subset of mesh cells with an appropriate size in postprocessing.Note that these arguments are also relevant for typical regular meshes where the size of the cells is variable.Although our examples include only energy spectra, the procedure is applicable for cross-spectra of two fields.Our main focus was on scale spectra E(k ℓ ), which are similar to the Fourier spectra, but the scale spectrum ϵ(ℓ) also deserves attention and might be similarly informative in practice.The procedure could be of interest for other analyses, such as computation of energy transfer across scales (see e.g., Aluie et al., 2018), but this raises an important question about the commutation of convolution and differential operations.Here we only outline the sources of potential difficulties.For staggered discretizations, the locations of discrete fields and the field derivatives are different.This does not create complications for regular quadrilateral meshes, where both sets of locations have the same patterns of neighbors, and the operations commute except at the boundaries.This is not the case even on uniform triangular or hexagonal meshes.For example, if a scalar field ϕ is defined at vertices of triangles (in a horizontal plane) as in FESOM, its gradient is located on triangles.If the convolution is applied first, it is operating in the space of values located at vertices.If one first computes gradients, and then applies the convolution, it operates in the space of values located on triangles.The discrete convolution operations in these two situations are different, and strict commutation is absent.The same problem will occur for other models (MPAS-Ocean, ICON-o), and also for other quantities such as velocity divergence.This problem will be absent for collocated discretizations.In addition, the implicit filters depend on the specified tolerance, the number of iterations, and preconditioners.The question on errors accompanying the exchange of the order of filtering and differencing operations on triangular meshes requires further studies.
The convergence of conjugate gradient solver slows down as k ℓ is approaching k min , and this slowdown becomes stronger when mesh resolution is increased.The biharmonic filter requires many more solver iterations than the harmonic one, and also needs a smaller tolerance.The convergence can be improved through the use of preconditioners.Guedot et al. (2015) used diagonal preconditioning.It provides an improvement for small k ℓ in our case too, but the number of iterations remains large (O(10,000)) for the biharmonic operator.In test implementation (in Matlab), the incomplete Cholesky preconditioner was efficient on regular meshes for the harmonic filter in essentially reducing the number of iterations (but not the execution time), however there are difficulties with using it on unstructured meshes or with the biharmonic filter.While the diagonal preconditioner and the conjugate gradient solver provide a functional implementation, it remains to be seen which preconditioners and solution methods will be able to improve the performance for general meshes.
Meshes with 1M wet vertices (the size of a typical quarter degree global mesh) can be treated in a serial way, requiring about 5 min for the construction of scale spectrum.A parallel implementation is needed for larger meshes.Both online and offline versions can be of interest.The online computations are straightforward for harmonic filters because they can rely on the existing solvers for the implicit sea surface height in FESOM and some other models.For biharmonic filters, the implementation in matrix form will require increasing halos to include neighbors of neighbors.A matrix-free form does not require the increased halos, but might complicate the design of preconditioners.Respective development is a subject of ongoing work.We did not try the higher-order filters (tri-harmonic, bi-biharmonic or higher) in this work.They might be needed to compute spectra that contain peaks.Technically, once discrete harmonic and biharmonic operators are available, any higher-order implicit filters can be realized as a sequence of more elementary inversions that use harmonic and biharmonic operators (Guedot et al., 2015).Practical implementation, convergence issues and the selection of the filter order need additional studies.
We experimented only with vertex placements for coarse-grained quantities, mostly because in this study we intended to explain and illustrate the concept.The cell (triangle) and edge placement will be addressed in future.In particular, for FESOM, which places horizontal velocities on triangles, the analysis of kinetic energy dissipation calls for computation on the original locations.Indeed, the projection of horizontal viscosity term from triangles to vertices may eliminate a significant part of small-scale dissipation, and computations on native locations will be made available in the future.The placement of coarse-grained velocities on triangles will lead to matrices with a twice larger dimension compared to the vertex case and a substantial increase in the required computational work.
Finally, as mentioned before, one of the advantages of the coarse-graining approach lies in the availability of spatial distributions of energy (or other quantities).These distributions allow one to learn not only which scales contain most of the energy, but also where this energy is mainly located.

Conclusions
We propose to use implicit filters to extend the concept of scale analysis based on coarse-graining to general unstructured meshes.The procedure is based on the discretizations of Laplacian operators and can provide filters of high-order if higher degrees of Laplacian are used.The discretization can be applied to scalar and vector data, both in flat and spherical geometry.
We demonstrated the effectiveness of the implicit coarse-graining approach through several examples, highlighting the similarity between scale spectra and standard Fourier spectra, as well as the ability to operate on diverse mesh types.
We acknowledge that numerous details regarding the extension of this approach to alternative data placements still need to be tested, which is the subject of ongoing research.
depicts these kernels as a function of z = kℓ = k/k ℓ .Note that the box kernel shows the sharpest behavior in the vicinity of zero, and the implicit second order filter is the most flat.The right panel shows the behavior of the dimensionless form-factor k ℓ ∂| Ĝℓ (k)| 2 /∂k ℓ = z∂ z | Ĝℓ (z)| 2 .The maxima are reached at z = ̅̅̅̅̅ 12

Figure 3 .
Figure 3.The side length of triangles (in km) calculated as 2A 1/ 2 c / 3 1/ 4 , where A c is the area of triangle c, on unstructured meshes U1 (left) and U2 (right).The two white circular regions correspond to land.The axes are in km.

Figure 4 .
Figure 4. Scale spectra computed with Laplacian smoothing (left) and biharmonic smoothing (right).Shown are the Fourier spectrum (red), the scale spectrum on a uniform triangular mesh obtained by bisecting quadrilateral cells (blue, mesh Q), on a regular equilateral mesh (solid black, mesh E) and on irregular unstructured meshes U1 and U2 (gray and green, respectively).The wavenumber is in cycle/km.The horizontal axis corresponds to K for the Fourier spectrum and to k ℓ = 1/ℓ for scale spectra.

Figure 5 .
Figure 5. Spectral peak detection by scale spectra based on biharmonic smoothing.The Fourier spectra are in red, the scale spectra are in blue.Peaks can be missed or masked if they are narrow and/or there are neighboring intervals with high spectral density.Section 2.3 gives estimates for the peak resolution in relation to the filter order.

Figure 6 .
Figure 6.Scale spectrum obtained by explicit coarse-graining with box filter (Equation12) (blue) compared to scale spectrum computed with Laplacian smoothing (solid black).The Fourier spectrum is in red, and the dashed line corresponds to the slope of 2.

Figure 7 .
Figure 7. Three variants of the same realization of ϕ and mesh in longitudelatitude coordinates (axes in degrees).

Figure 8 .
Figure8.Scale spectra computed for three meshes in Figure7using the biharmonic filter.The spectrum for the high-latitude variant is in green, and two other cases are under it.The red line corresponds to the Fourier spectrum on the parent regular quadrilateral mesh, and the dashed line gives the slope of 2.

Figure 9 .
Figure 9. Scale spectra of kinetic energy computed with the harmonic (left) and the biharmonic (right) filters.The latitude range in the mesh is approximately between 75 and 85°.The spectrum calculated with account for the metric terms is in blue, and the spectrum without the metric terms (in black) is only marginally visible at the smallest wavenumbers.The red line corresponds to the Fourier power spectrum of streamfunction multiplied with K 2 , and the dashed line indicates the slope of 2.