Metrics for Intercomparison of Remapping Algorithms (MIRA) protocol applied to Earth system models

Strongly coupled nonlinear phenomena such as those described by Earth system models (ESMs) are composed of multiple component models with independent mesh topologies and scalable numerical solvers. A common operation in ESMs is to remap or interpolate component solution fields defined on their computational mesh to another mesh with a different combinatorial structure and decomposition, e.g., from the atmosphere to the ocean, during the temporal integration of the coupled system. Several remapping schemes are currently in use or available for ESMs. However, a unified approach to compare the properties of these different schemes has not been attempted previously. We present a rigorous methodology for the evaluation and intercomparison of remapping methods through an independently implemented suite of metrics that measure the ability of a method to adhere to constraints such as grid independence, monotonicity, global conservation, and local extrema or feature preservation. A comprehensive set of numerical evaluations is conducted based on a progression of scalar fields from idealized and smooth to more general climate data with strong discontinuities and strict bounds. We examine four remapping algorithms with distinct design approaches, namely ESMF Regrid (Hill et al., 2004), TempestRemap (Ullrich and Taylor, 2015), generalized moving least squares (GMLS) (Trask and Kuberry, 2020) with post-processing filters, and WLS-ENOR (Li et al., 2020). By repeated iterative application of the high-order remapping methods to the test fields, we verify the accuracy of each scheme in terms of their observed convergence order for smooth data and determine the bounded error propagation using challenging, realistic field data on both uniform and regionally refined mesh cases. In addition to retaining high-order accuracy under idealized conditions, the methods also demonstrate robust remapping performance when dealing with non-smooth data. There is a failure to maintain monotonicity in the traditional L²-minimization approaches used in ESMF and TempestRemap, in contrast to stable recovery through nonlinear filters used in both meshless GMLS and hybrid mesh-based WLS-ENOR schemes. Local feature preservation analysis indicates that high-order methods perform better than low-order dissipative schemes for all test cases. The behavior of these remappers remains consistent when applied on regionally refined meshes, indicating mesh-invariant implementations. The MIRA intercomparison protocol proposed in this paper and the detailed comparison of the four algorithms demonstrate that the new schemes, namely GMLS and WLS-ENOR, are competitive compared to standard conservative minimization methods requiring computation of mesh intersections. The work presented in this paper provides a foundation that can be extended to include complex field definitions, realistic mesh topologies, and spectral element discretizations, thereby allowing for a more complete analysis of production-ready remapping packages.

1 Introduction

Coupled multimodel simulations often involve high degrees of computationally complex workflows, and achieving consistently accurate solutions is strongly dependent on the choice of spatiotemporal numerical algorithms used to resolve the interacting scales in physical models. Rigorous spatial coupling between components in such systems involves field transformations and communication of data across multiresolution grids while preserving key attributes of interest such as global integrals and local features, which is usually referred to as the process of remapping (also “regridding” or just “interpolation”) (Van Leer, 1979; Dukowicz and Kodis, 1987; Jones, 1999). Such remap procedures are critical in ensuring the stability and accuracy of scientific codes simulating multiphysics problems that typically occur in many different scientific domains. While there have been many high-order, stable interpolators proposed within different disciplines (Zienkiewicz and Zhu, 1992; Grandy, 1999; Jones, 1999; Smith et al., 2000; Garimella et al., 2007; de Boer et al., 2008; Slattery, 2016), to our knowledge, an effort to uniformly compare and measure the properties of these schemes has not yet been attempted. The current work is motivated by a need for an intercomparison of remapping schemes, which led us to standardize several numerical metrics to uniformly compare the properties of these algorithms that are routinely applied to problems in climate and weather system modeling.

Remapping algorithms, in general, compute the spatial interpolation or quasi-interpolation of field data that are defined on a source mesh (Ω_s) onto a target (Ω_t) component mesh. It is imperative to note that these mesh pairs, Ω_s and Ω_t, are generally of a different tessellation topology, spatial resolution, and regularity. Over the past several decades, there have been considerable efforts to create new conservative remapping algorithms to tackle various grid combinations and create consistently accurate schemes for solution transfer between component models in multiphysics problems in many fields besides climate. These efforts have resulted in several software libraries and packages developed for this purpose. Examples include schemes developed for fluid–structure interaction (FSI) or heat transfer (such as MpCCI, Joppich and Kürschner, 2006 and preCICE, Bungartz et al., 2016), moving mesh problems with arbitrary Lagrangian–Eulerian (ALE) methods (Dukowicz and Kodis, 1987; Dukowicz and Baumgardner, 2000), and general-purpose remap software such as MOAB (Tautges and Caceres, 2009; Mahadevan et al., 2020), PANG (Gander and Japhet, 2013), and Data Transfer Kit (DTK) (Slattery et al., 2013).

Remapping packages developed for ESMs, such as the Community Earth System Model (CESM) (Hurrell et al., 2013) and the Energy Exascale Earth System Model (E3SM) (E3SM Project, 2018), include SCRIP (Jones, 1999), the Earth System Modeling Framework (ESMF) Regridder (Hill et al., 2004), TempestRemap (Ullrich and Taylor, 2015; Ullrich et al., 2016), ncremap (Zender, 2008), and the Common Remapping Software (CoR) (Liu et al., 2018) libraries. These remapping packages generate files with mapping weights that are then consumed by couplers such as MCT (Larson et al., 2005), cpl7 (Craig et al., 2012), and OASIS3-MCT (Craig et al., 2017), which apply the weights during runtime. Support for generation of mapping weights, and application at runtime is growing (Mahadevan et al., 2020). Note that more recently, a comparative study of various production-ready remappers for ESMs using realistic meshes was performed by Valcke et al. (2022) using the specific numerical metrics that have been developed and presented in this paper. The success of this study in better understanding the numerical behavior of production regridders is promising.

Many of the production-ready remapping software implementations used in global climate simulations are typically based on first- and second-order conservative mesh-based schemes, with additional support for second-order nonconservative bilinear patch reconstructions (Zienkiewicz and Zhu, 1992; Rider, 2014) or third-order bi-cubic spline interpolations (Hanke et al., 2016) for scalar fields. The projection algorithms implemented in these libraries for the conservative maps are based on the computation of an overlay or a “supermesh” (Jiao and Heath, 2004a; Farrell et al., 2009; Farrell and Maddison, 2011), defined as ${\mathrm{\Omega }}_{S\cap T}={\mathrm{\Omega }}_{\mathrm{s}}\cap {\mathrm{\Omega }}_{\mathrm{t}}$, which is then used to compute consistent and conservative linear operators through L² minimization for transferring field information (Ullrich and Taylor, 2015; Ullrich et al., 2016).

Remapping also occurs in other parts of a climate model, such as within the geophysical fluid dynamics solver of a component. Remapping strategies for tracer transport advection schemes such as the Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme (Zerroukat et al., 2005), Conservative and Monotone Cascade Remap on Sphere for CubedSphere (CS) to regular latitude–longitude (RLL) grids (Lauritzen and Nair, 2008), geometrically exact conservative remapping approaches (Ullrich et al., 2009), and some other high-order mass-conserving schemes (Norman and Nair, 2008; Erath et al., 2013) have also been devised. While some of the previous work, including conservative semi-Lagrangian multi-tracer transport schemes (CSLAM) (Lauritzen et al., 2010) with second-order variants (Zerroukat et al., 2006), can be easily extended to arbitrary mesh topology combinations, many of the other approaches are tied to specific cases and not used in the inter-component coupling, which falls out of the scope of this study.

While traditional mesh-based schemes requiring computation of overlay–exchange grids have the advantage of being inherently conservative (Balaji et al., 2006), they tend to be computationally expensive (Jansen et al., 1992) and achieving high-order accuracy for scalar field data can sometimes pose difficulties. Numerically suboptimal approximations, such as nodal expansions or nearest neighbor values, can be used to compute interpolations quite efficiently if the field being remapped is smooth and does not mandate conservation prescriptions. Bilinear and bicubic interpolations can also prove to be alternatives here. However, consistently ensuring high-order accuracy under conservation constraints is often necessary for many real-world fields, such as the global heat flux that is exchanged between the atmosphere and ocean components in climate systems.

More recently, new mesh-based (Li et al., 2020) and fully meshless (Trask and Kuberry, 2020) remapping schemes that do not require overlap meshes have been developed for other problems in order to mitigate the computational complexity of computing exact intersections between unstructured meshes. These algorithms can provide high-order, conservative, and stable alternatives for climate modeling compared to the traditional linear maps generated from L²-minimization approaches.

As the number of available remapping algorithms grows, it becomes imperative to compare them under a unified framework to understand the properties of the schemes before applying them to real-world simulations. Additionally, while the computational cost is critically important for production runs, our intercomparison study specifically compares only the numerical performance of each algorithm under varying mesh topologies and field regularity, which are closely representative of those from a climate model such as E3SM. The presented protocol hence provides a systematic way to test and compare all existing and new remapping algorithms being developed for Earth system modeling.

Organization

The current study aims to better understand and document the key properties of the remapping schemes through the use of mesh, field, and scale resolution-independent metrics definitions. This intercomparison paper is organized as follows.

• A detailed literature survey of the current state-of-the-art remapping schemes used in climate modeling and various related coupled physics problems, along with the relevant numerical background for four specific high-order remapping algorithms and their implementations, is presented in Sect. 2.

• Definitions for remapping metrics that evaluate field accuracy, global conservation, strict global bounds control, and feature dispersion are presented in Sect. 3. We argue that these metrics are broadly useful for evaluating key properties that determine the accuracy and stability of the solution transfers between model components and can be applied to all remapping algorithms widely used in climate codes.

• A sample workflow for computing and comparing various numerical metrics in a unified and unbiased fashion is featured in Sect. 4. Some implementation details on the open-source Python-based infrastructure are also provided.

• Consolidated results from the intercomparison study applied to four competing remap algorithms on representative problems are shown in Sect. 5 for both uniformly refined and regionally refined global meshes.

• Potential future research directions to extend the analysis presented in this work are described in Sect. 6, followed by a summary of the intercomparison study in Sect. 7.

2 Remapping algorithms

In general, for coupled simulations, we need to transfer a field ${\mathit{\psi }}_{\mathrm{s}}\in {\mathbf{R}}^{n_{\mathrm{s}}}$ defined on the source grid, Ω_s, to a quantity ${\mathit{\psi }}_{\mathrm{t}}\in {\mathbf{R}}^{n_{\mathrm{t}}}$ on the target grid, Ω_t. A remap is hence defined as an operator $R:{\mathbf{R}}^{n_{\mathrm{s}}}\mapsto {\mathbf{R}}^{n_{\mathrm{t}}}$ that transfers ψ_s to a trial field ${\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{\mathrm{t}}\in {\mathbf{R}}^{n_{\mathrm{t}}}$. As noted above, we generally wish to preserve a property defined as an operator $P:{\mathbf{R}}^{n_{\mathrm{t}}}\to {\mathbf{R}}^m$ that describes discrete physical invariants and/or inequalities relevant to ψ_t. However, to keep the current work focused we restrict attention to definitions that are applicable only to scalar fields.

Finally, the accuracy or order of the remap typically depends on how the field ψ is reconstructed from a set of discrete values. On a mesh, this could be a Taylor series expansion around a cell centroid (Jones, 1999) or a general finite-element (FE) type of expansion using nodal basis functions.

2.1 Related work

One of the simplest data transfer methods is to use piecewise interpolation functionals. This approach is particularly convenient if ψ_s defined on Ω_s has an associated function space, which can be directly utilized for evaluating the interpolant approximation. Such consistent interpolation techniques utilizing the underlying basis functions for field descriptions give rise to standard second-order linear (on simplicial elements) or bilinear (with quadrilateral elements) interpolations (Hill et al., 2004). In finite-volume methods or for interpolating from station data (“scattered data” or “location streams”), the nearest-point interpolation is sometimes used, which is at best first-order accurate.

More generally, the remap methods designed for scattered data and cell-to-cell transfers are closely related to the reconstruction problem, i.e., to reconstruct a piecewise smooth function or its approximations at some discrete points on a mesh given some known quantities at discrete locations on the same mesh. The fundamental difference between reconstruction and remapping is that the former involves a mapping between a discrete and a continuous function space, while the latter involves a mapping between two discrete spaces defined on Ω_s and Ω_t. The demonstrated techniques used in high-order interpolators use spline quasi-interpolants (de Boor, 1990), bi-cubic splines (Hanke et al., 2016; Craig et al., 2017), the standard radial-basis function spaces (Flyer and Wright, 2007; Bungartz et al., 2016), the moving least squares (MLS) method (Lancaster and Salkauskas, 1981), and MLS variants such as the modified MLS (MMLS) (Joldes et al., 2015; Slattery, 2016), which originate from high-order reconstruction methods. Recent extensions to MLS for producing efficient high-order remap involve locally reconstructing the manifold geometry from a point set representation and then generating a compact stencil in the local coordinate chart (Liang and Zhao, 2013; Suchde and Kuhnert, 2019; Trask and Kuberry, 2020; Gross et al., 2020).

More recently, Li et al. (2020) proposed a high-order remap technique for piecewise smooth functions on surfaces, known as WLS-ENO remap (WLS-ENOR), which is based on the continuous moving frame (CMF) for smooth functions (Jiao and Wang, 2012) as well as an ENO-style technique (Liu and Jiao, 2016) for resolving discontinuities. For smooth functions, CMF differs from MLS (Lancaster and Salkauskas, 1981) in that it uses compact stencils over local coordinate systems from a 𝒞⁰ normal field instead of global stencils and 𝒞^∞ moving frames.

Typically, a remap method can and should be independent of the discretization methods used in physics models. In this context, the method may be nonconservative or conservative with respect to certain properties of the reconstruction. A major disadvantage of the quasi-interpolation techniques is the lack of inherent constraints for conserving energy during the remap process; see, e.g., de Boer et al. (2008) and Jiao and Heath (2004a). Note that the high-order quasi-interpolants (with degree of basis expansion p>1) tend to be significantly less dissipative than low-order methods ($p\in [\mathrm{0},\mathrm{1}]$); see, e.g., Slattery (2016). However, such remap methods are not guaranteed to produce conservative or monotone remapped solution fields.

A common approach to overcome this issue is to enforce conservation so that the integrals over the source and target meshes are equal. There are several examples of first- and second-order conservative remap schemes (Chesshire and Henshaw, 1994; Grandy, 1999; Gander and Japhet, 2013; Jones, 1999; Ullrich and Taylor, 2015) that rely on common-refinement-based L² projection (Jiao and Heath, 2004a) approaches. These schemes require computation of a function integral defined on the source mesh over some control volumes associated with a target node (or cell). The numerical integration is computed over the intersections of the elements (or cells) of the source and target meshes. These intersections form the common refinement (Jiao and Heath, 2004a), or supermesh (Farrell et al., 2009), whose computations require sophisticated computational geometry algorithms for efficiency and robustness (Jiao and Heath, 2004b, c). Although these first- and second-order schemes applied to ESMs are conservative, they may exhibit excessive numerical diffusion resulting in dissipation of “energy”, especially near field discontinuities or regions with large second derivatives.

Concurrently, when applying high-order remapping methods, discontinuities in the function defined on Ω_s can lead to overshoots or undershoots when evaluated on Ω_t. The discontinuities may be in the function values (also called 𝒞⁰ discontinuities), which tend to lead to 𝒪(1) oscillations (or “ringing”) that do not vanish under mesh refinement, analogous to Gibbs phenomena (Gottlieb and Shu, 1997; Fornberg and Flyer, 2007). Additionally, any discontinuities in the derivatives (𝒞¹ discontinuities) tend to cause milder oscillations that vanish under mesh refinement but nevertheless may accumulate in repeated remapping cycles. It is often critical to resolve or control these numerical artifacts as they introduce nonphysical variations that can influence the numerical stability of the coupled global nonlinear multiphysics system.

The deficiency in linear mapping approaches arises from the fact that they are only dependent on Ω_s and Ω_t but completely independent of the solution field that is being projected between the meshes. As a consequence of Godunov's theorem (Godunov, 1959) extended to linear remapping workflows with monotonicity constraints, there may be a restriction on the optimal achievable order of accuracy as shown by Van Leer (1979) while still preserving global solution bounds during the projection step. Hence, the properties of the linear maps can often be enhanced by using a procedure that is nonlinear (depending on the projected field variations) to recover property preservation for high gradient fields.

In this vein, the techniques for resolving Gibbs phenomena can be classified as filtering and mollification; see Jerri (2013). The filtering techniques often rely on post-processing, such as cropping and property redistribution, to ensure local conservation in a neighborhood. A recent variation of the mass borrowing approach (Royer, 1986) uses a limiter as a filter during the remapping process (Bradley et al., 2019) in order to impose local bounds preservation for linear map applications such that monotonicity can be recovered even when the underlying remapper does not provide it. This “Clip-And-Assured-Sum” (CAAS) post-processing filter can be useful to avoid spurious numerical oscillations due to resolution disparity or strong gradients in the underlying solution. Applying CAAS filters to quasi-interpolatory linear maps can hence produce strictly bounded reconstructions on Ω_t, while providing property preservation, as a viable option to achieve better remapped solutions in production climate simulations. In contrast, mollification adapts the kernels (i.e., basis functions) near discontinuities. In the past, discontinuity detecting and a posteriori stabilization procedures have been used with specific mesh discretizations (Blanchard and Loubere, 2016) to choose optimal orders of reconstruction adaptively in order to ensure better behavior for polygonal meshes. Other methods, such as WLS-ENOR, detect discontinuous regions and can effectively adapt the weighting schemes to resolve the Gibbs phenomena at the cost of additional computations at runtime.

Alternatively, rather than imposing a weakly nonlinear post-processing filter, using fully nonlinear remap schemes can be an option when the computational cost of the solution transfer is not the dominant factor in the simulation. Such nonlinear remap schemes typically use optimization-based remap (OBR) procedures (Carey et al., 2001; Bochev and Shashkov, 2005; Bochev et al., 2011) to minimize the net residual projection error using Lagrange multipliers. OBR follows a “divide-and-conquer” alternative (Bochev et al., 2014) to direct property preserving methods, which separates accuracy considerations from property preservation. In so doing, OBR helps to avoid interdependencies between mesh quality, accuracy, and property preservation that force many monotone reconstruction methods to make trade-offs between the latter two (Berger et al., 2005). Extensions of such schemes for remapping climate data are an unexplored topic and may provide avenues for future research.

Additionally, mimetic schemes that use compatible function spaces (Thuburn et al., 2009) depending on the fields being transferred between component models have proven to be extremely accurate for remapping scalar and vector fields on Arakawa C/D grids (Pletzer and Hayek, 2019). Potential extensions for arbitrary mesh topologies to yield compatible, conservative remaps for fluxes in climate components are possible (Ringler et al., 2010). However, to our knowledge, general theory and implementations for remapping on arbitrary meshes are currently unavailable, which may restrict the usage of such schemes for production climate simulations.

2.2 Weighted least squares approximations in remapping

A common theme across all of the remapping methods described in this paper is that they utilize some variants of the least squares approximations (also called quasi-interpolation) in their computational kernels. To illustrate the idea, let us consider a function $f\left(\mathit{u}\right):\mathrm{\Omega }\to \mathbb{R}$ at a given point ${\mathit{u}}_{\mathrm{0}}={\left[\mathrm{0},\mathrm{0}\right]}^T$, such as a quadrature point in a cell on the target mesh. Let us first assume a domain Ω⊂ℝ² for simplicity, where $\mathit{u}\equiv [u,v]$. Let f be 𝒞^p−1 continuous for some degree p≥0, and then f(u) can be approximated to p+1st-order accuracy about u₀ using a degree-p two-dimensional Taylor polynomial as

$$\begin{array}{}\text{(1)}& f\left(\mathit{u}\right)=\sum _{q=\mathrm{0}}^p\sum _{j,k\ge \mathrm{0}}^{j+k=q}{c}_{jk}{u}^j{v}^k+\mathcal{O}\left(h^{p+\mathrm{1}}\right),\end{array}$$

where $c_{jk}={\displaystyle \frac{\mathrm{1}}{j\mathrm{!}k\mathrm{!}}}{\displaystyle \frac{{\partial }^{j+k}}{\partial u^j\partial v^k}}f\left({\mathit{u}}_{\mathrm{0}}\right)$, and h is a measure of the radius of the local neighborhood. In cell-to-cell transfer, some integrals of f are typically known on the source mesh. Given m cells on the source mesh in a neighborhood of a point u₀ on the target mesh, let ϕ_i(u) be the test function (such as the Heaviside functions) associated with the ith cell. Let b_i denote these known integral values. We then obtain a system of m equations with $n=(p+\mathrm{1})(p+\mathrm{2})/\mathrm{2}$ unknown coefficients c_jk,

$$\begin{array}{}\text{(2)}& \underset{\mathrm{\Omega }}{\int }\sum _{q=\mathrm{0}}^p\sum _{j,k\ge \mathrm{0}}^{j+k=q}{c}_{jk}{u}_i^j{v}_i^k{\mathit{\varphi }}_i\left(\mathit{u}\right)\mathrm{d}\mathit{u}\approx b_i,\end{array}$$

for $\mathrm{1}\le i\le m$. Note that one can also use a bi-degree-p Taylor polynomial by letting $\mathrm{0}\le j\le p$ and $\mathrm{0}\le k\le p$ in Eq. (1), which would lead to $n=(p+\mathrm{1})(p+\mathrm{1})$ unknowns. The resulting approximate Taylor polynomial can then be used, for example, to evaluate (or quasi-interpolate) f at u₀ to p+1st order accuracy or even to higher order due to superconvergence. The same procedure can also be applied to construct a trivariate quasi-interpolation by replacing u and u^jv^k with x and x^jy^kz^ℓ, respectively. The m equations in Eq. (2) can be written in matrix form as Ax≈b, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ is known as a generalized Vandermonde matrix, x∈ℝⁿ is composed of c_jk in Eq. (1), and b∈ℝ^m is composed of the known integrals b_i on the source mesh. In general, the generalized Vandermonde system (2) is rectangular. It can be solved by minimizing a weighted norm of the residual vector $\mathit{r}=\mathit{b}-\mathbf{A}\mathit{x}$, i.e.,

$$\begin{array}{}\text{(3)}& \underset{\mathit{x}}{min}\parallel \mathit{r}{\parallel }_{\mathbf{W}}\equiv \underset{\mathit{x}}{min}\parallel \mathbf{W}(\mathbf{A}\mathit{x}-\mathit{b}){\parallel }_{\mathrm{2}},\end{array}$$

where W is a diagonal matrix containing the weight for each row in A. This formulation is known as weighted least squares or WLS for short (Golub and Van Loan, 2013). When m=n and A and W are nonsingular, W does not affect the solution. More generally, when m≠n, different W leads to different solutions by changing the relative importance of certain sample points. Different remapping algorithms and reconstruction schemes use specific weighting strategies to achieve optimal solutions to the weighted least squares problem in Eq. (3). Note that the generalized Vandermonde matrix A may be ill-conditioned as the polynomial degree p used for the reconstruction grows. A preferred approach to address this potential ill-conditioning is to use a rank-revealing QR factorization (Golub and Van Loan, 2013; Li et al., 2020).

2.3 Algorithms for Earth system models

Among the standard conservative remapping and high-order reconstruction strategies introduced in Sect. 2.1 for climate problems, we selected four specific algorithms to explore in detail. The motivation for choosing these algorithmic implementations was driven by their high usability, including the use in current ESMs, and the rigorous underlying numerics that can be verified and validated consistently for a large suite of test problems. We categorize these algorithms below into three broad groups based on whether the algorithms require overlay meshes and whether mesh data structures are utilized to compute the solution reconstruction.

• I.

  Overlay-mesh-based remappers

  We consider two specific implementations that provide conservative remapping capability for production ESMs.

  • a.

    ESMF: the Earth System Modeling Framework's Regrid function providing bilinear and conservative maps

  • b.

    TempestRemap: conservative, consistent, and monotone remapper with higher-order L² projection support

  Both ESMF Regrid and TempestRemap provide the remapping capability for scalar fields defined on Ω_s based on a weighted residual formulation given in Eq. (3). With a source field ψ_s defined on Ω_s, the formulation computes ψ_t on Ω_t by solving the following problem:

  $$\begin{array}{}\text{(4)}& \underset{\mathrm{\Omega }}{\int }{\mathit{\psi }}_{\mathrm{t}}{\mathit{\varphi }}_i\mathrm{d}V=\underset{\mathrm{\Omega }}{\int }{\mathit{\psi }}_{\mathrm{s}}{\mathit{\varphi }}_i\mathrm{d}V,\end{array}$$

  where the ϕ_i represents suitable weight functions. In a common-refinement-based transfer (Jiao and Heath, 2004a), the ϕ_i represents the basis functions ψ_t,i on Ω_t, which leads to a Galerkin projection method. Such a Galerkin method can achieve conservation globally by solving the weighted residual minimization in Eq. (4). More details on the specific methodologies used for reconstruction in these implementations are provided in later sections.

  Note that the SCRIP library (Jones, 1999) also provides both first- and second-order conservative remapping schemes, but it was not included in the current study since its capabilities are similar to that of ESMF and TempestRemap implementations, and those two are used in current production versions of CESM and E3SM.

• II.

  Meshless remappers

  Reconstruction methods that do not directly utilize the topological information about the underlying mesh layout are meshless methods. In our study, we will consider the generalized moving least squares (GMLS) method, with global conservation, monotonicity constraints, and local bounds preservation provided by CAAS as a post-processing filter.

  Future studies could also include the comparison of MMLS (Slattery, 2016) and RBF (Bungartz et al., 2016) reconstructions within this framework for remapping climate data, since those methods have demonstrated some success in fluid–structure interaction and nuclear engineering problems.

• III.

  Non-overlay mesh-based hybrid remappers

  The final category includes the mesh-based remappers that do not require computation of an intersection mesh between Ω_s and Ω_t. The weighted least squares essentially non-oscillatory remap method (WLS-ENOR) utilizes a discontinuity-capturing, high-order technique for piecewise smooth functions to produce optimal field transformations between meshes by minimizing the residual in Eq. (3).

  Other opportunities for comparison in this category include using reconstructed climate data with the conservative bilinear algorithm and patch reconstructions in ESMF or bicubic interpolations available in Yet Another Coupler (YAC) (Hanke et al., 2016).

These algorithms and their implementations span a range of remapping techniques, including those currently used in production runs and those that can potentially be used in the future given the availability of open-source software. Further details regarding the numerics and the software tools for each of the schemes are provided in the following subsections.

2.3.1 ESMF

The Earth System Modeling Framework (ESMF, https://earthsystemmodeling.org/, last access: 7 August 2022) (Hill et al., 2004) consists of a suite of software modeling tools that support building Earth system models. Among other features, ESMF exposes several key functionalities to transfer data between component models in weather and climate applications. It offers a variety of data structures for transferring field data between components and libraries for regridding, time advancement, and other common modeling functions. The advanced regridding algorithms provided by ESMF implement standard bilinear interpolation, higher-order patch recovery, first- and second-order conservative projections, and several variations of nonconservative nearest neighbor interpolants. ESMF can produce maps that are “offline” in the sense that they can be precomputed, stored, and then applied as a linear operators to arbitrary fields defined on Ω_s to compute the projection on Ω_t. Fully online remapping is also possible with ESMF.

In the current intercomparison study, we utilized the command-line applications installed with ESMF (version 8.1) to generate and apply interpolation weights from the command line using NetCDF files. Using this ESMF_RegridWeightGen application, interpolation weights can be generated for various grid combinations and applied to source field data to compute projections between component models. More specifically, we select only the first- and second-order conservative projection methods invoked with “conserve” and “conserve2nd” command-line options, respectively, which are implemented in ESMF for scalar fields and fluxes. These options are routinely used in production climate models and hence provide valuable baselines on the current state of remapping algorithms in our comparison study. For more details regarding the numerics and implementation of the algorithms in ESMF, we refer readers to Collins et al. (2005), Jones (1999), and Kritsikis et al. (2017).

2.3.2 TempestRemap

TempestRemap (https://github.com/ClimateGlobalChange/tempestremap, last access: 7 August 2022) is a software library to generate conservative, consistent, and monotone remapping weights of arbitrarily high-order accuracy (with degree p>1) between meshes on the sphere. Similar to the ESMF tools, TempestRemap generates offline maps that are consumed by ESMs. TempestRemap supports the computation of remapping weights between any combination of finite-volume (FV) and/or spectral element (SE) (both continuous and discontinuous spaces) descriptions with conservative prescriptions. For purposes of the current paper, we describe only the high-order FV–FV algorithms below.

In TempestRemap, the procedure used to generate remapping weights for FV discretizations consists of two primary operations (Ullrich and Taylor, 2015; Ullrich et al., 2016). First, local polynomial reconstructions are defined over the source mesh with some adjustments made for spherical geometry (Jalali and Gooch, 2013). The coefficients of these local reconstructions are computed according to a weighted pseudo-inverse method (Skamarock and Menchaca, 2010; Skamarock and Gassmann, 2011). These polynomials are integrated over the target mesh by using the overlap, or supermesh (Farrell et al., 2009), which in general can be defined as ${\mathrm{\Omega }}_{\mathrm{s}\cap \mathrm{t}}={\mathrm{\Omega }}_{\mathrm{s}}\cap {\mathrm{\Omega }}_{\mathrm{t}}$. Note that when computing high-order linear maps, if a sufficiently large patch on the source map is not used during the reconstruction process, the condition number of the generalized Vandermonde matrix A in Eq. (2) can be very high, leading to numerical roundoff errors and degradation in the overall accuracy of the remap operator to first-order accuracy in the neighborhood. Further details on the analysis are provided in Ullrich and Taylor (2015).

A common approach for the integration operator has been to construct a potential function with divergence equal to the scalar field being integrated and then use the divergence theorem to transform integrals over mesh faces into line integrals around the face (Dukowicz and Kodis, 1987; Jones, 1999; Ullrich et al., 2009). While this technique has been used to define conservative remapping operators, a drawback is that it requires an analytical expression for the potential function, which can be difficult in general. An alternative method was proposed in Erath et al. (2013), which, although avoiding some of the difficulties of the line integral approach, lacks consistency. In TempestRemap, a quadrature-based integration operator is used to generate an initial set of remapping coefficients, which is guaranteed to be consistent (exactly remaps the unit field) but may not be conservative. To obtain conservation, the coefficients of the matrix are projected linearly into the space of maps that are both conservative and consistent.

2.3.3 GMLS

The generalized moving least squares (GMLS) method extends the moving least squares (MLS) technique from approximation of point values to approximation of arbitrary linear functionals (Nayroles et al., 1992; Breitkopf et al., 2000; Wendland, 2004; Mirzaei et al., 2012). Similar to MLS, a distance-based weighting kernel is used to favor the source data sites closer to the point of query or reconstruction. GMLS features a broad choice of sampling functionals and target operators, and hence the term “generalized” is in its name. The application of nonlinear sampling functionals and target operators for GMLS is possible (Gross et al., 2020) but is certainly not common. In fact, most applications of GMLS use linear sampling functionals and target operators (gradient, curl, divergence, integral average along an edge or over a cell, etc.), for which there is theory on approximation and well-posedness. When the sampling functionals and the target operator are simply pointwise evaluations, GMLS reduces to the traditional MLS method.

High-order accurate approximations cannot be achieved with traditional MLS schemes using function spaces described by Raviart–Thomas (ℋ(div)) and Nedelec basis (ℋ(curl)) or even cell-averaged quantities that are common in FV codes. However, through careful selection of the sampling functionals and target operators, consistent approximations of fields embedded in these various nonstandard spaces are possible using the GMLS approach. This flexibility of the GMLS method greatly increases the available types of input data that can be handled to produce high-order reconstructions of fields between Ω_s and Ω_t. When a GMLS problem is solved, a set of coefficients corresponding to elements of a basis (traditionally polynomial) is computed. The target operation is then applied to each member of the basis used for approximation, after which a dot product is made between the originally computed coefficients and the evaluation of the target operation acting on the basis. While GMLS permits a wide choice of target operators, we choose a target operation which is the cell average of each member of the approximating basis on the target grid.

Traditionally, GMLS uses a basis that is defined as a function of the spatial dimension from which a point cloud is sampled. However, in this work, reconstruction of functions sampled on a manifold permits generating a compact stencil in a local coordinate chart, which is one dimension smaller (Liang and Zhao, 2013; Suchde and Kuhnert, 2019; Trask and Kuberry, 2020; Gross et al., 2020). The savings in net computational floating-point operations (flops) is a factor of p³ in R² compared to a traditional basis in R³, where p is the order of the basis.

As is true for many other regression schemes, GMLS is not inherently conservative to machine precision, but rather it is “conservative” to discretization precision. In other words, the degree to which it violates conservation is discretization-dependent and generally vanishes with refinement. However, such weak conservation notions may not be deemed satisfactory for climate modeling applications for which exact global conservation has a history of being demanded and valued. In such cases, GMLS remap should be followed by a post-process filter to restore global conservation to machine precision. In this paper, we use either the GMLS or GMLS-CAAS notations to indicate whether the CAAS routine has been used as a post-processing filter after each remap step in order to restore global conservation or global bounds, along with an attempt to improve local property preservation.

Similar to the overlay-mesh-based, ESMF, and TempestRemap offline remappers described in Sect. 2.3.1 and 2.3.2, respectively, the solution of the GMLS problem produces a stencil or a linear mapping that can be computed a priori to a simulation. This enables storage of the relevant parts of the remap operator as a sparse matrix. Therefore, the use of GMLS without post-process filtering can be thought of as an offline remapper, and GMLS-CAAS can combine the offline and online processes to achieve more favorable properties. Note that the CAAS filter is only one of several choices available for the online post-processing filter. Alternative strategies such as the nonlinear OBR (Bochev et al., 2014) for feature detection and minimization of local dissipation may be possible in conjunction with the GMLS workflow. These enhanced GMLS remapping strategies may be considered in future studies.

There are several key motivations for using GMLS to perform field remapping for ESMs. These include mesh topology independence, as well as flexibility in the choice of the sampling functional and the target operator, thereby enabling remap of nontraditional degrees of freedom that may be defined on the vertices, edges, or faces of Ω_s. Additionally, the embarrassingly parallel nature of the dense linear systems that require inversions during reconstruction yields a high flops-to-communication ratio. This feature makes GMLS more suited for next-generation platforms. Given the computationally intensive nature of the GMLS algorithm, it is best implemented in libraries focused on parallel performance portability, computational efficiency, and significant compiler optimizations. The Compadre Toolkit (https://github.com/SNLComputation/compadre, last access: 7 August 2022) (Kuberry et al., 2019) is built on Kokkos (Edwards et al., 2014) and Kokkos kernels such that a single version of the code is written to be performant on both multicore CPU and hybrid GPU architectures. Compadre is built to assemble and solve large batches of GMLS problems in parallel, thereby leveraging batched QR solvers with pivoting for the parallel solution of many small dense linear systems.

2.3.4 WLS-ENOR

WLS-ENOR, or weighted-least-squares-based essentially non-oscillatory remap (Li et al., 2020), is based on the WLS-ENO reconstruction technique proposed in Liu and Jiao (2016). Originally developed for solving hyperbolic partial differential equations (PDEs), WLS-ENO detects discontinuities and then reduces or eliminates the potential Gibbs phenomena in the solutions of the PDEs by adapting the weighting schemes in WLS based on the local features in the solution. WLS-ENOR adapted WLS-ENO to remap data between meshes, and in the process, it enhanced the treatment of discontinuities to resolve not only the severe oscillations at 𝒞⁰ discontinuities (i.e., jumps in function values), but also the accumulated effect of mild oscillations due to 𝒞¹ discontinuities (i.e., jumps in the derivatives).

Unlike the preceding remapping techniques, WLS-ENOR is a non-overlay mesh-based technique in that it uses the mesh for computing numerical integration, but it does not require an overlay mesh (although it has the option to use the overlay mesh if available). More specifically, WLS-ENOR utilizes adaptive quadrature rules with p refinement in smooth regions and h refinement near discontinuities in its numerical integration to achieve high-order accuracy (with degree p>1) and (nearly local) conservation. More specifically, WLS-ENOR is composed of three major components. The first component is a WLS-based algorithm for smooth functions, as we described in Sect. 2.2. In this context, the weighting scheme in WLS-ENOR is based on a positive-definite radial function due to Buhmann (Buhmann, 2001). As shown in Li et al. (2020), this weighting scheme encourages statistical error cancelations and enables better superconvergence (i.e., higher than (p+1)st-order convergence) with even-degree p for node-to-node transfer of smooth functions. In this intercomparison work, we use an extension of the work in Li et al. (2020) to cell-centered data. Mathematically, this extension essentially uses the step functions (also called Heaviside functions) over the cells on the source mesh as the test functions in Eq. (2) compared to the Dirac delta function at the nodes on the source mesh as test functions in Li et al. (2020). Note that in this work, we apply WLS on a sphere Ω=S, which is topologically a 2D object embedded in ℝ³. One could construct a WLS in ℝ³ directly as in Slattery et al. (2013), which would lead to a large generalized Vandermonde system. Alternatively, we can construct a local uv (surface tangent) coordinate frame at each point x₀∈S. Specifically, let m₀ denote an approximate normal at x₀, which can be the exact normal to S or a first-order estimation. Let t₁ and t₂ form an orthonormal basis of an approximate tangent plane orthogonal to m₀. The local uv coordinate frame is then centered at x₀ with axes t₁ and t₂. We can then transform the sample points about x₀ to use this local uv coordinate frame and apply the WLS to construct a bivariate quasi-interpolation. In terms of implementation, WLS-ENOR constructs a matrix-based transfer operator for smooth functions, which maps the cell-averaged values from the source mesh to the target mesh.

The second component in WLS-ENOR is the detection and resolution of discontinuities. In particular, WLS-ENOR can detect 𝒞⁰ and 𝒞¹ discontinuities of the function on the source mesh and then transfer the discontinuity tags from the source mesh to the target mesh. The detector in WLS-ENOR is based on an asymptotic analysis of the Gibbs phenomena near 𝒞⁰ and 𝒞¹ discontinuities as described in Li et al. (2020). The original detector in Li et al. (2020), however, was for node-to-node transfer. For cell-to-cell transfer in this work, we simply apply the detector to the dual mesh and treat cell-averaged values as approximations of cell-center values. This approximation is second-order accurate, which is sufficient in detecting discontinuities. After identifying the discontinuities, WLS-ENOR uses a solution-based weighting scheme that effectively leads to (nearly) one-sided quasi-interpolation in discontinuous regions. This solution-based weighting scheme causes WLS-ENOR to use a different set of basis functions to overcome the Gibbs phenomena, and hence it can be classified as a mollification technique (see Sect. 2.1). In terms of implementation, the discontinuity detector on the source mesh is implemented as a matrix-based operator; after the discontinuities are identified, a new local solution-based transfer operator is constructed for each cell on the target mesh near a discontinuous cell on the source mesh. We refer readers to Li et al. (2020) for details.

The third component in WLS-ENOR is an adaptive quadrature technique, which is enabled when the target mesh is significantly coarser than the source mesh. In this case, simply sampling the function at the quadrature points of a target cell may miss some important local features on the source mesh, especially near discontinuities. To overcome this loss of information, one could use the overlay mesh as in TempestRemap. This approach may be ideal in terms of accuracy and conservation, but it introduces complications when the elements have curved edges. Although WLS-ENOR implementation supports this option, in this comparative work, we use a non-overlay-based version of WLS-ENOR that utilizes h and p refinement of the quadrature rules in discontinuous and smooth regions, respectively. In particular, near the detected discontinuities, WLS-ENOR subdivides the cells (i.e., h-refinement) on the target mesh to match the local resolution on the source mesh and then utilizes the quadrature points of the subdivided cells. For smooth regions, we found it more efficient to use the quadrature points of higher-degree quadrature rules (i.e., p refinement) than subdividing the cells. This adaptive quadrature technique not only overcomes the loss of accuracy but also enables WLS-ENOR to recover global conservation in a fashion that is nearly local to discontinuous regions. To this end, if there is a gain or loss in terms of global conservation, we distribute this global error proportionally to cells that have gained or lost mass locally, correspondingly. This conservation–recovery procedure reduces the pollution of the conservation errors from discontinuous regions to smooth regions. We estimate the local conservation errors using the subdivided cells near discontinuities for accuracy and robustness; for smooth regions, we use a simple comparison of the local extreme values in the local neighborhood for smooth regions for efficiency.

The current implementation of the WLS-ENOR algorithm uses MATLAB, with which the core components were converted into C using “MATLAB Coder” (version 4.2). An open-source C++ implementation is currently underway and will be released in the future for both node-to-node and cell-to-cell field transfers.

3 Metrics to evaluate remapping algorithms

Solution remapping on unstructured meshes is a complex process, and it is critical to satisfy several key properties to ensure that the transferred field data between components do not introduce unbounded and nonphysical error modes. In order to compare different remapping schemes in an unbiased framework, we introduce five primary categories under which the comparison metrics can be grouped.

1. Sensitivity is algorithmic invariance to underlying component mesh topology.

2. Consistency is the ability to retain the order of convergence of the underlying discretization in a given norm.

3. Conservation entails ensuring global integral (and) local conservation of critical quantities.

4. Monotonicity involves the preservation of global and local solution bounds during remap and solution transfer between components.

5. Dissipation is the minimization of local solution dispersion on repeated back-and-forth remap transfers between Ω_s and Ω_t.

Given a continuous field ψ, we use D^s and D^t to denote the reference discretization of ψ on the source and target grid, respectively. These are generated by directly sampling the analytical fields and by spherical harmonic expansion (see Sect. 4) for the realistic fields. The regridding operator from the source to target mesh is denoted by R; i.e., the regridded field is denoted by RD^s[ψ]. The global integral operator is denoted by I^s and I^t on the source and target grid, respectively. Typically these operators take the form

$$\begin{array}{}\text{(5)}& {\displaystyle }I\left[\mathit{x}\right]=\sum _{\text{all }j}{x}_j{J}_j,\end{array}$$

where J_j denotes the area or weight of the jth degree of freedom (i.e., the volume of a finite volume).

3.1 Grid sensitivity

A crucial factor for the success and broader usability of a general remapping algorithm in ESMs is the ability to produce mesh-independent numerical behavior that is robust for any pair of structured or unstructured meshes. In other words, remapping algorithms need to be general and without approximations targeted at specific topological elements.

In the current work, we utilized three different meshes of varying resolutions. Specifically, we present analysis performed to compare remapping schemes using the cubed-sphere (CS), quasi-uniform Voronoi (MPAS), regular latitude–longitude (RLL) meshes on both quasi-uniform and regionally refined meshes (RRMs). Some sample meshes used in the study are shown in Fig. 1.

Figure 1 A depiction of the five meshes studied in validation, unit testing, and intercomparison of the regridding schemes.

3.2 Standard accuracy measures

Accuracy in the remapped solution will be assessed with standard error metrics defined as follows.

$$\begin{array}{}\text{(6)}& {\displaystyle }{∥E∥}_{L_{\mathrm{1}}}& {\displaystyle }\equiv {\displaystyle \frac{I^{\mathrm{t}}\left[{\mathit{e}}_K\right]}{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}\right]}},\text{(7)}& {\displaystyle }{∥E∥}_{L_{\mathrm{2}}}& {\displaystyle }\equiv \sqrt{{\displaystyle \frac{I^{\mathrm{t}}\left[{\mathit{e}}_K^{\mathrm{2}}\right]}{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]{\mathrm{|}}^{\mathrm{2}}\right]}}},\text{(8)}& {\displaystyle }{∥E∥}_{L_{\mathrm{\infty }}}& {\displaystyle }\equiv {\displaystyle \frac{max\left[{\mathit{e}}_K\right]}{max\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}}},\end{array}$$

where ${\mathit{e}}_K=\mathrm{|}{\mathbf{RD}}^{\mathrm{s}}\left[{\mathit{\psi }}_K\right]-{\mathbf{D}}^{\mathrm{t}}\left[{\mathit{\psi }}_K\right]\mathrm{|}$ is the local remapping error in element K∈Ω_t, E represents the discrete error vector in the solution field relative to the reference field sampled on Ω_t, and I^t is the weighted integral using Eq. (5) on Ω_t. In some general sense, the error measure ${∥E∥}_{L_{\mathrm{1}}}$ identifies errors in large-scale features, ${∥E∥}_{L_{\mathrm{2}}}$ identifies errors in small-scale features, and ${∥E∥}_{L_{\mathrm{\infty }}}$ identifies the largest pointwise error. These are special cases of error measures ${∥E∥}_{L_p}$ in L_p norm, in which as p tends from 1 to ∞, the norms capture the large-scale mean absolute errors and the largest pointwise errors at the two extremes, respectively, for p=1 and p=∞. Related to accuracy, consistency is assessed by applying these metrics to uniformly smooth fields with no C⁰ or C¹ discontinuities and verifying the asymptotic theoretical rate of convergence under uniform mesh refinement conditions. The asymptotic convergence order of a given remapping method with a degree-p reconstruction is in general expected to be 𝒪(h^p+1) in ${∥E∥}_{L_{\mathrm{2}}}$ and 𝒪(h^p) in ${∥E∥}_{H_{\mathrm{1}}}$, where h is the characteristic spatial length of the mesh that can be computed using a simple definition such as $h=\frac{\mathrm{1}}{N_{{\mathrm{\Omega }}_{\mathrm{t}}}\left(\mathrm{el}\right)}$, with $N_{{\mathrm{\Omega }}_{\mathrm{t}}}\left(\mathrm{el}\right)$ representing the number of elements in the unstructured mesh Ω_t. Note that for nonuniform meshes such as RRMs, it is important to take the root mean square (rms) element Jacobian value for the parameter h. Additionally, while the accuracy of the remap projection is strongly dependent on both Ω_s and Ω_t, in studies presented here, we only utilize uniform refinements in both meshes (without including cross-resolution calculations resulting from a fine Ω_s and coarse Ω_t) to compute the overall convergence order for projecting smooth fields.

Note that in order to eliminate potential aliasing errors, the normalization factors D^t[ψ] used in the denominator for definitions of ${∥E∥}_{L_{\mathrm{1}}}$, ${∥E∥}_{L_{\mathrm{2}}}$, and ${∥E∥}_{L_{\mathrm{\infty }}}$ are computed based on the exact sampling (element-averaged for FV discretization) of the data (reference solution) on Ω_t, and not using the projection of the field RD^s[ψ]∈Ω_t.

3.3 Gradient preservation measures

Preservation of the solution gradients in addition to other critical properties, such as local conservation in the remapping procedure, requires C¹ continuity in the D^s[ψ]. Let ∇D^s and ∇RD^t be the gradients of the scalar fields on Ω_s and Ω_t, respectively.

Then, in order to measure accuracy of the solution and its gradient, we introduce two specific global metrics: ${\left|E\right|}_{H_{\mathrm{1}}}$ semi-norm and the ${∥E∥}_{H_{\mathrm{1}}}$ norm.

$$\begin{array}{}\text{(9)}& {\displaystyle }{\left|E\right|}_{H_{\mathrm{1}}}^{\mathrm{2}}& {\displaystyle }\equiv {\displaystyle \frac{I^{\mathrm{t}}\left[\mathrm{\nabla }{\mathit{e}}_K\mathrm{\nabla }{\mathit{e}}_K\right]}{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]{\mathrm{|}}^{\mathrm{2}}\right]}},\text{(10)}& {\displaystyle }{∥E∥}_{H_{\mathrm{1}}}^{\mathrm{2}}& {\displaystyle }\equiv {\displaystyle \frac{I^{\mathrm{t}}\left[{\mathit{e}}_K^{\mathrm{2}}+\mathrm{\nabla }{\mathit{e}}_K\mathrm{\nabla }{\mathit{e}}_K\right]}{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]{\mathrm{|}}^{\mathrm{2}}\right]}}={∥E∥}_{L_{\mathrm{2}}}^{\mathrm{2}}+{\left|E\right|}_{H_{\mathrm{1}}}^{\mathrm{2}},\end{array}$$

where ${\mathit{e}}_K=\mathrm{|}{\mathbf{RD}}^{\mathrm{s}}\left[{\mathit{\psi }}_K\right]-{\mathbf{D}}^{\mathrm{t}}\left[{\mathit{\psi }}_K\right]\mathrm{|}$ is the local remapping error in element K, $\mathrm{\nabla }{\mathit{e}}_K=\mathrm{|}\mathrm{\nabla }{\mathbf{RD}}^{\mathrm{s}}\left[{\mathit{\psi }}_K\right]-\mathrm{\nabla }{\mathbf{D}}^{\mathrm{t}}\left[{\mathit{\psi }}_K\right]\mathrm{|}$ is the corresponding gradient of the error, and I^t is the weighted integral using Eq. (5) on Ω_t.

In this study, ∇q for some scalar q, defined as a cell mean value, is generated by finding the convex hull of cells surrounding the current cell and computing the gradient per Barth and Jespersen (1989). We compute these gradients as part of the metrics evaluation after each remapping sequence.

3.4 Global conservation

Global conservation is trivially assessed by evaluating the change in the global integral of the scalar field value on the source mesh and the projected field on the target mesh. We use the following metric to quantify global conservation.

$$\begin{array}{}\text{(11)}& {\displaystyle }{L}_g\equiv {\displaystyle \frac{I^{\mathrm{t}}\left[{\mathbf{RD}}^{\mathrm{s}}\left[\mathit{\psi }\right]\right]-I^{\mathrm{s}}\left[{\mathbf{D}}^{\mathrm{s}}\left[\mathit{\psi }\right]\right]}{I^{\mathrm{s}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{s}}\left[\mathit{\psi }\right]\mathrm{|}\right]}}\end{array}$$

However, we note that this definition for L_g is only meaningful when the target domain fully envelops the source domain (which may have gaps or holes in more general cases). In the case of climate modeling, an admissible example for using Eq. (11) would be for remapping heat and moisture fluxes from the land surface to the overlying atmosphere.

3.5 Global extrema preservation

Global extrema preservation can be assessed via the standard G_min and G_max error metrics (Ullrich and Taylor, 2015):

$$\begin{array}{}\text{(12)}& {\displaystyle }\left|G_{\min }\right|& {\displaystyle }\equiv min\left\{\mathrm{0},{\displaystyle \frac{min\left({\mathbf{RD}}^{\mathrm{s}}\left[\mathit{\psi }\right]\right)-min\left({\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\right)}{max\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}-min\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}}}\right\},\text{(13)}& {\displaystyle }\left|G_{\max }\right|& {\displaystyle }\equiv max\left\{\mathrm{0},{\displaystyle \frac{max\left({\mathbf{RD}}^{\mathrm{s}}\left[\mathit{\psi }\right]\right)-max\left({\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\right)}{max\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}-min\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}}}\right\}.\end{array}$$

The error measures G_min and G_max identify undershoots and overshoots, respectively, by taking on nonzero values ($\left|G_{\min }\right|\le \mathrm{0}$ and $\left|G_{\max }\right|\ge \mathrm{0}$) when there is a departure away from the reference global extreme values. In other words, a nonzero value of the metric indicates changes in global extrema, indicating the presence of a Gibbs phenomenon or instabilities introduced due to the interpolation. Hence, the global extrema metric is particularly useful as it provides indications about the monotonicity-preserving properties of the remapping schemes.

3.6 Local extrema preservation

Local extrema preservation can be assessed using a localized difference; i.e., to what degree does the remapped grid cell value fall within the range of surrounding grid cells sampled on the target grid? This consideration motivates us to define a localized difference in extrema:

$$\begin{array}{ll}{\displaystyle }{\mathrm{\Delta }}_{\min ,j}& {\displaystyle }\equiv min\left\{\mathrm{0},{\left({\mathbf{RD}}^{\mathrm{s}}\left[\mathit{\psi }\right]\right)}_j\right.\\ \text{(14)}& {\displaystyle }& {\displaystyle }\left.-\underset{\text{1-ring patch}}{min}{\left({\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\right)}_j\right\},\forall j\in {\mathrm{\Omega }}_{\mathrm{t}},{\displaystyle }{\mathrm{\Delta }}_{\max ,j}& {\displaystyle }\equiv max\left\{\mathrm{0},{\left({\mathbf{RD}}^{\mathrm{s}}\left[\mathit{\psi }\right]\right)}_j\right.\\ \text{(15)}& {\displaystyle }& {\displaystyle }\left.-\underset{\text{1-ring patch}}{max}{\left({\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\right)}_j\right\},\forall j\in {\mathrm{\Omega }}_{\mathrm{t}},\end{array}$$

where the minimum and maximum are taken over all grid cells on the target mesh that surround the current target cell j. These values can then be reduced to a single value in the usual manner by using an appropriate norm definition for both Δ_min,j and Δ_max,j:

$$\begin{array}{}\text{(16)}& \begin{array}{l}{L}_{min,\mathrm{1}}\equiv \frac{I^{\mathrm{t}}\left[\mathrm{|}{\mathrm{\Delta }}_{\min }\mathrm{|}\right]}{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\right[\mathit{\psi }\left]\mathrm{|}\right]},\\ L_{\min ,\mathrm{2}}\equiv \frac{\sqrt{I^{\mathrm{t}}\left[{\mathrm{\Delta }}_{\min }^{\mathrm{2}}\right]}}{\sqrt{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]{\mathrm{|}}^{\mathrm{2}}\right]}},\\ L_{min,\mathrm{\infty }}\equiv \frac{{max}_j\mathrm{|}{\mathrm{\Delta }}_{\min ,j}\mathrm{|}}{max\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}-min\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}}.\end{array},\begin{array}{l}{L}_{max,\mathrm{1}}\equiv \frac{I^{\mathrm{t}}\left[\mathrm{|}{\mathrm{\Delta }}_{\max }\mathrm{|}\right]}{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\right[\mathit{\psi }\left]\mathrm{|}\right]},\\ L_{max,\mathrm{2}}\equiv \frac{\sqrt{I^{\mathrm{t}}\left[{\mathrm{\Delta }}_{\max }^{\mathrm{2}}\right]}}{\sqrt{I^{\mathrm{t}}\left[\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]{\mathrm{|}}^{\mathrm{2}}\right]}},\\ L_{max,\mathrm{\infty }}\equiv \frac{{max}_j\mathrm{|}{\mathrm{\Delta }}_{\max ,j}\mathrm{|}}{max\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}-min\mathrm{|}{\mathbf{D}}^{\mathrm{t}}\left[\mathit{\psi }\right]\mathrm{|}}.\end{array}\end{array}$$

Note that the definition of the localized differences shown in Eqs. (14) and (15) utilizes a local neighborhood to determine the deviation from reference extrema values. This is sufficient to capture resolution of sharp gradients in the remapped fields under mesh refinement for element-averaged data. However, the metric contains 𝒪(h) dependence on the mesh resolution and can be applied to 𝒞⁰ or smoother fields, but not when 𝒞⁰ discontinuities are present.

4 Metrics for Intercomparison of Remapping Algorithms (MIRA) workflow

For all remapping algorithms evaluated in this comparison study, we conduct iterative two-way (Ω_s→Ω_t and Ω_t→Ω_s) remapping of an initial source field with FV discretization on Ω_s. We do so in order to characterize the stability of a scheme and expose any dissipation effects, which would not be possible to ascertain when comparing single, unidirectional field transfers. With this workflow, we seek to quantify the consistency, stability, and convergence of each participating algorithm as measured with the metrics defined in Sect. 3. Here, N_r indicates the number of iterative applications of the linear map to compute field transformation on ${\mathrm{\Omega }}_{\mathrm{s}}\to {\mathrm{\Omega }}_{\mathrm{t}}\to {\mathrm{\Omega }}_{\mathrm{s}}$, which will be referred to as remap iterations. While production ESM solvers do not utilize repeated remap transfers at every time step, our approach to use iterative remaps can provide valuable insight into the dissipative effects and long-term temporal behavior of fully coupled climate simulations, in addition to determining the stability of the remap operator without a need for explicit spectral analysis.

4.1 Open-source MIRA implementation

The workflow necessary to evaluate a given remapping method comprises five consecutive steps described below.

1. Generate a series of meshes of different topologies and resolutions. We use the cubed-sphere (CS), quasi-uniform Voronoi (MPAS), regular latitude–longitude (RLL), and regionally refined (RMM) grids of the CS and/or MPAS types of varying resolutions to devise the test cases. See Fig. 1 for an illustration of the meshes. The mesh data are stored in universal NetCDF4 format containing an array of vertex point locations and a cell connectivity map to describe the topology.

2. Given a collection of meshes as in step 1 above, a Python module called MeshPreprocessDriver is then used to generate and store the adjacency maps and unstructured cell area integrals with high-order Gauss quadrature rules. The convex hull map for each cell is also precomputed and stored during this step in order to speed up the evaluation of remapped field gradient metrics.

3. A second Python module called FieldGenDriver then takes each of the pre-processed mesh files and evaluates scalar fields by sampling from either an analytical function on the sphere or a set of prescribed spherical harmonic (SPH) coefficients, which is described in Sect. 4.2. In this step, a cell average is computed by local quadrature within each mesh element to a given order of accuracy by appropriately choosing the order of the quadrature to resolve the SPH expansion order. This operation is performed on all the input meshes to generate the reference “ground truth” realization of a given field, which is used to accurately compute the metrics defined in Sect. 3.

   We emphasize that any existing mesh (such as Yin–Yang Kageyama and Sato, 2004, or cubic–octahedral – Gaussian – reduced grids) can be used in this workflow, instead of the ones generated in step 1, as FieldGenDriver only relies on the existence of element connectivities and adjacency maps to be available for computing cell integrals.

4. All remapping algorithms evaluated in this study use the mesh data, depending on the scheme, and initial reference solutions on Ω_s to execute the test suite over one or many N_r iterations. The expected outputs from each of the algorithms for the test problems devised are the discrete solution vectors ${\mathit{\psi }}_{\mathrm{s}}^i\in {\mathrm{\Omega }}_{\mathrm{s}}$ and ${\mathit{\psi }}_{\mathrm{t}}^i\in {\mathrm{\Omega }}_{\mathrm{t}}$, where $i\in [\mathrm{1},\mathrm{\dots }{N}_{\mathrm{r}}]$. In the current study, unless otherwise specified, N_r=1000.

5. The final Python module in the metrics suite, MetricsDriver, can then be invoked on each of the remapped output data (typically stored in NetCDF files) to consistently compute all the remapping metrics defined in Sect. 3. The computed remapping metrics are then stored as comma-separated values (CSVs) for further analysis and intercomparison studies.

The schematic shown in Fig. 2 provides further details on this workflow. Note that in order to evaluate and compare a new remapping implementation such as SCRIP (Jones, 1999), CoR (Liu et al., 2018), or YAC (Hanke et al., 2016), only the fourth and fifth steps in the workflow have to be executed, since the pre-processed input meshes and the sampled reference data have been made publicly available (Mahadevan et al., 2021).

Figure 2 The MIRA workflow for generating the remapping metrics for the intercomparison study.

4.2 Scalar test variables on the sphere

In the remapping intercomparison study, we consider five scalar test variables defined on the sphere as reference solutions fields. These fields are chosen such that different aspects of the remapper can be evaluated uniformly. Details about the analytical and real-world fields as well as the sampling methodology used in the Python implementation are provided below.

4.2.1 Idealized analytical fields

Idealized fields used in this study mirror the approach of (Lauritzen and Nair, 2008) and (Ullrich and Taylor, 2015). Namely, we employ three idealized test cases of varying complexity to understand the error measures produced by remapping. The two analytical fields studied are depicted in Fig. 3.

The first analytical field (AnalyticalFun1) is a combination of spherical harmonics functions with frequency wave similar to order 3, given by

$$\begin{array}{}\text{(17)}& \mathit{\psi }=Y_{\mathrm{3}}^{\mathrm{2}}+Y_{\mathrm{3}}^{\mathrm{3}},\end{array}$$

where $Y_m^l$ represents the real spherical harmonic functions evaluated through the SHTools package for degree m and polynomial order l.

Following (Jones, 1999), and (Lauritzen and Nair, 2008), the second field (AnalyticalFun2) is a relatively smooth function resembling spherical harmonics of order 2 and azimuthal wavenumber 2, given by

$$\begin{array}{}\text{(18)}& \mathit{\psi }=\mathrm{2}+{\cos }^{\mathrm{2}}\mathit{\theta }\cos \left(\mathrm{2}\mathit{\lambda }\right),\left(Y_{\mathrm{2}}^{\mathrm{2}}\right).\end{array}$$

These fields are used to test performance for a smooth, well-resolved field and a slightly high-frequency, weakly resolved field with rapidly changing gradients. Given that the analytical expressions for these fields are trivial to evaluate, we can compute the exact numerical errors introduced by the remapping schemes when projecting the fields from Ω_s to Ω_t.

Note that the FieldGenDriver module can take arbitrary closed-form functions and evaluate them on the sphere by using high-order quadrature order rules to sample and compute element-averaged data. This design allows the flexibility to test slightly more complex analytical vortex fields (Ullrich et al., 2009) or any three-dimensional real-valued function projected on the sphere with coordinate transformations (Townsend et al., 2016).

Figure 3 Contour plots of analytical fields used in this study: AnalyticalFun1 (a) and AnalyticalFun1 (b).

4.2.2 Real data fields

We also test the performance of each remap technique by regridding real data fields obtained from freely available composite satellite observations. The fields chosen are total precipitable water (TPW), cloud fraction (CFR), and global topography (TOPO). From a representative dataset, we compute a spherical harmonic (SPH) decomposition in order to determine an analytical approximation of spectral content. We employ these particular fields because we can control their characteristics as functions, i.e., global bounds for topography, positive definiteness for precipitable water, and continuity for cloud fraction. As such, these fields present distinct challenges to the remapping methods. For convenience, we use the SHTools (Wieczorech and Meschede, 2018) package through its Python interface (PySHTools), which facilitates the computation of spectra based on spherical harmonic bases and reconstruction of fields thereof.

Given the 1D averaged spatial amplitude spectrum for each set of composite satellite data as shown in Fig. 4 with a corresponding linear fit, we then produce controlled randomized realizations for each field on any unstructured mesh and resolution, including regionally refined grids. The randomization is applied to the coefficients of the expansion at each degree, extracted from the linear fit functions in Fig. 4, and is entirely reproducible given an integer seed to a pseudo-random number generator (set to 384 in our study). The code for this is available in the FieldGenDriver (Guerra et al., 2021). The original data range from 0.25 to 0.1^∘ of resolution and SPH reconstructions are smoothly varying up to the number of coefficients employed. Note that while the order of SPH reconstruction can be specified arbitrarily high (between 1 and 512 modes) to get a better resolution of the fields, the computation of element-averaged sampling representative of FV discretization needs to utilize a sufficiently high-degree quadrature rule such that the SPH expansions are exactly integrated on the sphere.

Figure 4 1D averaged spatial amplitude spectrum for real data fields based on composite satellite observations. Fit coefficients over two branches (dotted blue and green lines) correspond to a function $f=a{x}^b+c$, where x is the logarithm of the spherical harmonic degree.

Figure 5 Randomized reconstruction of the real-world fields used in this study.

Total precipitable water (TPW).

Global composite data for TPW are taken from the MIMIC-TPW2 project (http://tropic.ssec.wisc.edu/real-time/mtpw2/product.php?color_type=tpw_nrl_colors&prod=global2&timespan=24hrs, last access: 7 August 2022) (Wimmers and Velden, 2011). Morphological compositing is applied to microwave sensor data and numerically advected. This field has an absolute minimum for all realizations of 0.0 mm with maxima typically in the 70.0 to 80.0 mm range. A representative, randomized reconstruction of TPW is shown in Fig. 5 (top-left).

Cloud fraction (CFR).

Global composite data for CFR are taken from the NASA AQUA/MODIS data archive (https://neo.gsfc.nasa.gov/view.php?datasetId=MYDAL2_M_CLD_FR, last access: 8 August 2022) (Platnick et al., 2020). This field uses absolute global limits for all realizations of 0.0 to 1.0. Thus, by imposing these bounds after each reconstruction, 𝒞⁰ discontinuities are introduced. These manifest as spatially flat regions in the data where maximum cloud cover is noted. A randomized reconstruction of CFR is shown in Fig. 5 (top right).

Global topography (TOPO).

Global topography data are taken from the ETOPO1 Global Relief Model (https://www.ngdc.noaa.gov/mgg/global/global.html, last access: 7 August 2022) (NOAA National Geophysical Data Center, 2009; Amante and Eakins, 2009). This field has a global minimum and maximum for all realizations of −10994.0 and 8848.0 m but is otherwise smoothly varying. A randomized reconstruction of TOPO is shown in Fig. 5 (bottom).

The raw satellite data snapshot of TPW, CFR, and TOPO fields that are used in the current study have been made available separately (Guerra and Mahadevan, 2021) in order to make the workflow reproducible.

5 Results and discussion

Comparing different remapping algorithms under a unified infrastructure for test problems and metrics collection is a nontrivial task. The metrics defined in this study and the implementation of the various field samplings on arbitrary unstructured meshes have provided large output datasets to analyze the key properties of the remappers under consideration.

Specifically, for uniformly refined experiments, a series of different mesh types (CS, MPAS, RLL) $N_{\mathrm{type}}^{\mathrm{uni}}=\mathrm{3}$ with five different refinement levels ($N_{\mathrm{ref}}^{\mathrm{uni}}=\mathrm{5}$) was generated. Similarly, for the regionally refined experiments, two different mesh types ($N_{\mathrm{type}}^{\mathrm{rrm}}=\mathrm{2}$) using CS and MPAS with three refinement levels ($N_{\mathrm{ref}}^{\mathrm{rrm}}=\mathrm{3}$) around the continental US were used. Table 1 provides the details on the number of elements and nodes in the various meshes utilized in the current study.

Table 1 Details on the number of elements (N_Ω(el)) and nodes (N_Ω(nd)) for different mesh types and resolutions used in the study: (a) uniformly refined meshes and (b) regionally refined meshes.

Using the field definitions (N_fields=5) introduced in Sect. 4.2, consisting of two smooth analytical fields and three real fields expanded with spherical harmonics, sampling was then evaluated on all input meshes ($N_{\mathrm{type}}^{\mathrm{uni}}{N}_{\mathrm{ref}}^{\mathrm{uni}}+N_{\mathrm{type}}^{\mathrm{rrm}}{N}_{\mathrm{ref}}^{\mathrm{rrm}}=\mathrm{15}+\mathrm{6}=\mathrm{21}$) to serve as reference solutions. With these input meshes (Mahadevan et al., 2021), the remapping algorithms were applied following the workflow in Fig. 2 for combinations of CS-MPAS (both uniform and RRM), MPAS-RLL, and RLL-CS meshes of varying resolutions to evaluate the metrics data.

The volume of consolidated output metrics data is enormous from this experiment, since 1000 remapping iterations were performed on $N=N_{\mathrm{type}}^{\mathrm{uni}}{C}_{\mathrm{2}}[N_{\mathrm{ref}}^{\mathrm{uni}}{]}^{\mathrm{2}}{N}_{\mathrm{fields}}=\mathrm{375}$ global, uniformly refined mesh cases and $N=N_{\mathrm{type}}^{\mathrm{rrm}}{C}_{\mathrm{2}}[N_{\mathrm{ref}}^{\mathrm{rrm}}{]}^{\mathrm{2}}{N}_{\mathrm{fields}}=\mathrm{45}$ RRM cases, with each of the four remapping algorithms using various degrees of reconstructions p for FV–FV field transfers. Measuring more than 15 different remapping metrics for each of these cases has provided extensively detailed results to compare the algorithmic implementations in an unbiased fashion. Hence, unless explicitly noted, only significantly unique results are presented and discussed below in the following subsections. We direct readers to the IPython notebooks available in Mahadevan et al. (2021), which can be used to generate the comparison for any combination of mesh types, source resolution, target resolution, field variable, and remapping schemes. Note that all the metrics data collected during the analysis are also stored in the same repository for reproducibility.

Detailed results from the intercomparison study and discussion on the implication of each metric to the remapping scheme are presented next.

5.1 Consistency

The consistency of the high-order remap algorithm implementations can be verified by remapping smooth functions and calculating the spatial convergence order of the resultant approximations on the target mesh after repeated remaps. Using the sampled analytical functions described by AnalyticalFun1 and AnalyticalFun2, a verification study was conducted using the standard error norms and gradient error metrics data for the various schemes. For smooth solution profiles, the theoretically expected convergence rates for a consistent remapping method are 𝒪(h^p+1) in ${∥E∥}_{L_{\mathrm{1}}}$, ${∥E∥}_{L_{\mathrm{2}}}$, and ${∥E∥}_{L_{\mathrm{\infty }}}$ global error norms and 𝒪(h^p) in ${∥E∥}_{H_{\mathrm{1}}}$ and ${\left|E\right|}_{H_{\mathrm{1}}}$ global gradient error norms. The convergence rates measured in the current studies use the definition of h provided in Sect. 3.2 by using the N_Ω(el) data described in Table 1 for different refined meshes.

In general, from these studies we observed that convergence rates for both low- and high-degree approximations up to p=2 of the analytically smooth fields show good agreement with the theoretically expected accuracy convergence rates. However, for high-order mesh-based remaps, achieving a convergence rate higher than 𝒪(h³) appears to be a limitation with the methods involved in this study, while the hybrid and meshless schemes show consistent recovery of the approximation orders for all degrees tested. The following subsections provide detailed results and discussions for each of the remapping algorithms.

5.1.1 ESMF

The conservative schemes implemented in the ESMF package (Hill et al., 2004) have been thoroughly verified and are routinely used to generate the linear maps for solution transfer between components in E3SM. The second-order conservative projection algorithm that was originally introduced for ALE computations (Dukowicz and Kodis, 1987), and later applied to spherical meshes (Jones, 1999), has been implemented in ESMF with an appropriate linear gradient reconstruction (Kritsikis et al., 2017). We measure the convergence rates for both the first-order (conserve) and second-order (conserve2nd) conservative schemes and present the results observed in Table 2. The ESMF first-order conservative scheme yields expected rates of 𝒪(h) asymptotically. However, the second-order scheme (conserve2nd) shows degraded convergence rates as confirmed by the global error and gradient norms. This convergence result is unexpected and contrary to the analysis of the second-order, piecewise linear finite-volume reconstruction procedure presented by Kritsikis et al. (2017), which is implemented in ESMF.

Table 2 ESMF: convergence rates for the AnalyticalFun1 field on all mesh types.

Note that we have presented the computed convergence rates from the analysis as is, and qualitatively speaking, the values in Table 2 that are approximately equal to 1.0 indicate 𝒪(h) and approximately equal to 0.0 indicate 𝒪(1) behavior.

In order to better understand the relative accuracy of the first- and second-order conservative remapping implementations in ESMF, we performed a comparative analysis for all grid combinations using standard global error norms. These results are shown in Fig. 6, where the ${∥E∥}_{L_{\mathrm{2}}}$ and ${∥E∥}_{H_{\mathrm{1}}}$ error profiles for both the conservative schemes are compared as a function of remap iterations (N_r) for the r4 resolution meshes. The computed error metrics from both the schemes indicate that, even though the convergence rates are similar, the conserve2nd option in ESMF produces a significantly better approximation for the remapped field, irrespective of the mesh resolution or field characteristics in the tested samples. It is also evident that the ESMF algorithms are more accurate for the RLL-CS mesh types in all global error norms, indicating sensitivity to the underlying element types in the mesh. The significantly better error approximations of the conserve2nd conservative scheme emphasize that it should be preferred over the conserve scheme when possible, irrespective of the underlying mesh types involved in remap.

Figure 6 ESMF comparison of conserve and conserve2nd options for all mesh type combinations (CS-MPAS, MPAS-RLL, RLL-CS) of resolution r4 and the TOPO field using log(${∥E∥}_{L_{\mathrm{2}}}$) (a) and log(${∥E∥}_{H_{\mathrm{1}}}$) (b) global error metrics as a function of N_r.

5.1.2 TempestRemap

The conservative high-order linear maps computed by TempestRemap, as shown in Table 3, produce higher-order expected theoretical convergence rates in comparison to ESMF for the smooth analytical solution fields. The convergence results presented here have been generated with a bi-degree-p basis reconstruction using a rectangular truncation strategy that computes source patches based on edge-adjacency graphs.

Note that even for smooth solutions, the convergence rates observed for p>2 degrade during the conservative reconstruction for FV–FV maps. This degradation is evident in the rate reduction in all the global error norms for p>2, which indicates that a larger patch size may be necessary to accurately and consistently recover the smooth field after remap from Ω_s to Ω_t. The rate reduction is especially noticeable in the presence of pole singularities like that of RLL meshes, where the ${∥E∥}_{L_{\mathrm{\infty }}}$ shows further reduction in convergence rates. The failure to achieve 𝒪(h⁴) or higher-order accuracy is an artifact of the implementation in TempestRemap, which can be improved with further analysis and is not a limitation of the underlying numerical method.

Table 3 TempestRemap: convergence rates for the AnalyticalFun1 field on all mesh types.

TempestRemap produces conservative solution projections between mesh combinations that can be third-order with p=2 for smooth fields. However, even if local element-wise conservation can be guaranteed in overlay-based high-order L²-minimization schemes, monotonic reconstructions may not be strictly possible without additional effort. This behavior is because Godunov's theorem (Godunov, 1959) precludes the existence of optimal high-order linear maps that are also monotone. Due to this restriction, and since Gibbs phenomena (Jerri, 2013) are ubiquitous with high-order maps when steep gradients are encountered, global bounds can be preserved in TempestRemap by employing nonlinear filtering algorithms such as CAAS during the online solution transfer in ESMs.

5.1.3 GMLS

The remapping scheme based on the meshless generalized moving least squares (GMLS) method demonstrates the flexibility to deliver higher-order convergence for scalar field data. The convergence rates computed for the AnalyticalFun1 field are shown in Table 4 for various polynomial degrees.

Table 4 GMLS and GMLS-CAAS: convergence rates for the AnalyticalFun1 field on all mesh types. Note that the superconvergence of gradients for p=4 appears to be special for AnalyticalFun1; for AnalyticalFun2, GMLS converged at about fourth order, as theoretically expected, in ${∥E∥}_{H_{\mathrm{1}}}$ and ${\left|E\right|}_{H_{\mathrm{1}}}$.

The convergence rates for the nominal GMLS scheme and the GMLS-CAAS remapping method with a post-processing step in Table 4 show that, in general, high-order accuracy can be achieved. However, using the nominal GMLS scheme for climate modeling problems can result in nonconservative and potentially oscillatory reconstructions for fields with strong gradients. Hence, GMLS-CAAS algorithm is especially advantageous to enable global and local bounds preservation. Note that the augmented GMLS-CAAS algorithm suffers from convergence degradation for higher polynomial degree values and is limited to 𝒪(h³) for this smooth field data in ${∥E∥}_{L_{\mathrm{2}}}$. There is limited theoretical proof for convergence rates of the CAAS filter in the literature (Bradley et al., 2019) for arbitrary problems. However, in two dimensions, one can consider a field with 1D connected bands of extrema that demonstrates a maximum rate of convergence of 𝒪(h³) in ${∥E∥}_{L_{\mathrm{2}}}$ norm. Algorithms like CAAS function by clipping newly formed local extrema resulting from higher-order reconstructions and computing a redistribution of the mass deficit accordingly. We hypothesize that the observed convergence degradation is primarily a result of these clipping and redistribution steps.

Even though the GMLS-CAAS remaps show lower convergence rates, it provides the benefit of making the scheme globally conservative and monotone. The CAAS algorithm requires runtime modification of the projected fields to ensure global and local bounds, and the nonlinear solution-dependent filter can eliminate Gibbs oscillations, providing better stability during remap operation.

5.1.4 WLS-ENOR

The convergence rates for various polynomial degrees of reconstruction p are tabulated in Table 5. The convergence analysis for the WLS-ENOR scheme shows that even for high polynomial degrees, theoretically expected rates of 𝒪(h^p+1) in ${∥E∥}_{L_{\mathrm{2}}}$ and 𝒪(h^p) in ${∥E∥}_{H_{\mathrm{1}}}$ are observed for the smooth analytical fields. For example, for p=4, we can confirm that the asymptotic convergence rate for ${∥E∥}_{L_{\mathrm{2}}}$ is 𝒪(h⁵) and for ${∥E∥}_{H_{\mathrm{1}}}$, ${\left|E\right|}_{H_{\mathrm{1}}}$ is 𝒪(h⁴).

Table 5 WLS-ENOR: convergence rates for the AnalyticalFun1 field on all mesh types.

We note that the WLS-ENOR algorithm is equipped with an internal nonlinear filtering (or more precisely, mollification) mechanism to detect sharp gradients and discontinuities in order to adaptively choose the weights during the high-order reconstruction process locally. In contrast to the GMLS-CAAS high-order meshless scheme with a post-processing filter that results in a convergence order degradation, the WLS-ENOR scheme remains consistently accurate for smooth functions up to p=4 in our experiments.

5.1.5 Real fields: convergence rate comparisons

While high-order convergence rates are achievable for smooth field profiles of AnalyticalFun1 and AnalyticalFun2, maintaining theoretical rates of convergence for “real-world” field data is not guaranteed because such data may lack the regularity necessary to achieve the best possible rates. Due to the presence of strong 𝒞⁰ and 𝒞¹ discontinuities in the real scalar fields, which are representative of flux fields exchanged between components in coupled climate solvers, a remap comparison on the SPH-sampled TPW field data was computed, as tabulated in Table 6. The convergence rate comparisons using different remap schemes show that the computed rates for higher-order polynomial reconstructions fall severely short of theoretical rates as expected relative to the rates shown previously for AnalyticalFun1.

However, it can be observed that higher than first-order schemes tend to show better behavior on all meshes and fields tested. Additionally, the linear maps computed with ESMF underperform the other schemes, as evident from the gradient error norms ${∥E∥}_{H_{\mathrm{1}}}$ and ${\left|E\right|}_{H_{\mathrm{1}}}$. The meshless schemes are competitive in terms of accuracy bounds and provide a viable alternative approach for usage in production climate simulations, in which L²-minimization schemes and low-order bilinear maps have traditionally been routinely used.

Table 6 Comparing global error measures for the TPW field with different remap schemes on all mesh types of finest (r4) Ω_s and Ω_t resolutions.

Note that in this particular comparison, we selected the highest polynomial degree p that was tested for each remapping method, even though for production climate simulations, p=1 or p=2 may typically yield a more numerically stable solution due to Gibbs phenomena. In such operational circumstances, high-order methods for TempestRemap and GMLS should be augmented with nonlinear post-processing filters, such as CAAS.

5.2 Global conservation

All mesh-based L² projection scheme implementations chosen in this study, namely ESMF and TempestRemap, are globally conservative by the nature of the underlying numerics (Ullrich and Taylor, 2015). Since these implementations compute consistent reconstructions using overlay meshes, any deficit with respect to exact conservation can be verified by performing column sum operations on the discrete linear map matrix representing the projection or mass-matrix operator.

The global field integral error indicating a conservation deficit for the TPW field variable is shown in Fig. 7 for different mesh resolutions of the CS-MPAS mesh type. In these plots, the y axes represent the logarithmic error of the global integral value between the solution remapped onto Ω_t relative to the reference solution sampled on Ω_t. Hence, lower values indicate better satisfaction of the global conservation metric.

The expected behavior of preserving global integrals in the remapped solutions, as illustrated in Fig. 7, has been verified for both low- and higher-order maps, irrespective of the mesh topology or the field variables being transferred across meshes. Notably, the accumulation of roundoff error in the global field integral as a function of remap iterations can quickly add to the global conservation error when using linear maps. This effect is seen in the field projections computed with both TempestRemap and ESMF through two sparse matrix–vector (SpMV) products in each remap iteration. As expected based on computational complexity of the SpMV operation, the low-order conservative ESMF remaps do accumulate errors more slowly than high-order TempestRemap projections, especially at fine mesh resolutions.

Figure 7 Global conservation error metric for the TPW field on CS-MPAS meshes with different resolutions in semi-log scale.

The WLS-ENOR scheme achieves excellent global conservation by applying adaptive quadrature rules to integrate the functions to (nearly) machine precision and redistributing the conservation errors locally near discontinuous regions. The unmodified GMLS scheme is nonconservative without any post-processing and hence is not presented here. However, using the nonlinear CAAS filtering algorithm with the GMLS remapping scheme provides global conservation to user-specified tolerances.

This conservation metric describes whether the remapped solution preserves the global integral over the domain, which is often a desirable constraint for most scalar variables in order to ensure that the total mass and energy in the closed simulation system remain constant. However, for some scalar fields such as sea surface temperature (SST) or TPW, such a conservation constraint may not be strictly mandated, and this enables the use of even the nonconservative GMLS scheme, which has been demonstrated to achieve excellent accuracy.

5.3 Monotonicity – global extrema preservation

Departures away from global extrema, i.e., overshoots and undershoots, provide a useful metric to assess the stability of the remap operator under investigation. A perfectly stable method has a zero value for this metric. The results for the monotonicity metric as a function of the remap iterations are shown in Figs. 8 and 9 for $\left|G_{\max }\right|$ of the CFR variable on MPAS-RLL and $\left|G_{\min }\right|$ of the TOPO variable on CS-MPAS meshes, respectively. These field variables have been specifically chosen to showcase the effects of the remapping algorithms on the well-defined upper and lower global field bounds that cannot be violated.

The results clearly demonstrate that the mesh-based remapping schemes such as those in ESMF and TempestRemap implementations do not adhere to the global extrema (both maxima and minima) after the first linear remap operator application. The magnitude of this departure is resolution-dependent, but the behavior is consistent in all cases tested. The use of low-order, dissipative linear maps from ESMF shows a slow drift away from global maxima, which are nearly recovered after several remaps (around 700 in Fig. 8). However, the global minima increase continuously within the bounds of the reference solution. Hence, the conserve2nd option in ESMF damps the global maxima and amplifies the global minima. While not currently pursued in the paper, further analysis of the spectral properties of the low-order forward and reverse maps can shed light on the dampening properties of these ESMF linear maps.

In contrast, the high-order TempestRemap solutions show a drift away from the bounds in every case. The mildly damped increase in the metric as a function of remap operator application signals the presence of spurious high-frequency modes in the linear map, which is more prevalent as the degree of reconstruction increases. Obtaining high-order remapped solutions in addition to preserving monotonicity in these schemes will require post-processing filters such as CAAS, which can be used to improve the behavior in TempestRemap for such problems of interest.

Compared to the traditional remapping schemes, the meshless GMLS-CAAS and hybrid WLS-ENOR remap schemes maintain the global bounds for most test cases and fields tested in this study. Augmenting the non-monotone GMLS scheme with the CAAS bounds preservation algorithm shows excellent, stable behavior for all cases. WLS-ENOR shows relatively high compliance with actual bounds because its non-oscillatory weights near discontinuities essentially preserve the monotonicity and convexity of the solutions. However, when the meshes are sufficiently fine, the functions may appear smooth relative to the grid resolution, and in such scenarios, WLS-ENOR would not apply the adaptive weights and the global extremes may oscillate randomly at a magnitude comparable to the local discretization error, as seen in Fig. 8b. Note that the hidden traces of $\left|G_{\min }\right|$ and $\left|G_{\max }\right|$ for GMLS-CAAS and WLS-ENOR in Figs. 9 and 8, respectively, are at the expected levels of zero (monotone reconstructions). This is evident for the GMLS-CAAS and WLS-ENOR methods in many of the experiments.

Figure 8 $\left|G_{\max }\right|$ metric for the CFR field on MPAS-RLL with different target mesh resolutions.

Figure 9 $\left|G_{\min }\right|$ metric for the TOPO field on CS-MPAS with different source mesh resolutions.

5.4 Locality – local extrema preservation

Similar to the global bounds metrics that measure solution monotonicity, the local bounds metrics compute the norm of error resulting from the comparison of the reference sampled field data against the projected field in a one-ring local neighborhood. This metric provides insight into the preservation of local field features on repeated remap operator applications. The computed locality metrics for maxima and minima for various mesh resolutions and real-world fields are shown in Figs. 10 and 11.

These results showcase the fact that the high-order TempestRemap implementations perform quite well with minimal feature loss for all fields on all meshes. In contrast, the low-order ESMF linear maps have high dissipation growth that is dependent on both the resolution and remap iterations. For high-resolution meshes, the amplification of the dissipation at every iteration can even be undamped and grows linearly with the ESMF conserve2nd option. This result implies that as the number of remap iterations increase, more information on local features is lost during every transfer for these cases.

For the WLS-ENOR scheme, which uses a hybrid strategy to detect local features and adaptive quadrature rules to resolve a high discrepancy in mesh resolutions, its preservation of local bounds exhibits a strong dependence on the spatial resolution of Ω_s and Ω_t. However, unlike ESMF, the slope of metric indicating departure away from local bounds remains near zero as remap iterations are increased, signifying that the reconstructions remain stable on repeated applications without any further feature losses, probably due to the use of adaptive quadrature rules.

In the meshless remapping category, the GMLS schemes utilize the CAAS algorithm to clip and conserve fields in a local subset of the neighborhood for each evaluation point. For GMLS-CAAS, local bounds must be given to CAAS, and it is important to consider how they are determined. As GMLS is a meshless solution technique, local bounds for each site of reconstruction are determined by computing the maximum and minimum of the values in the neighborhood (determined meshlessly using a K–d tree) used for reconstruction from Ω_s. With a coarse Ω_s resolution, these local bounds computed from K–d tree queries are no longer accurate estimators to enforce tight bounds in a one-ring neighborhood of Ω_t, which is the metric measured in this study for local bounds preservation. The discrepancy in bounds for such scenarios could be minimized by using a mesh intersection of the cells in the source and target neighborhoods, but no attempt was made to do this for GMLS-CAAS in keeping with the spirit of its application as a purely meshless technique.

Figure 10 $L_{max,\mathrm{2}}$ metric for the TPW field on CS-MPAS meshes with different resolutions.

Lauritzen and Nair (2008) noted that when remapping using their monotone and conservative CaRS algorithm, the higher-order reconstructions do not significantly improve accuracy when Ω_s is finer than Ω_t. However, the reverse case shows significant benefit with high-order polynomial reconstructions. These conclusions are similar to those observed in the current intercomparison experiments shown in Figs. 10 and 11. Note that in Fig. 10c, for equal-resolution CS-MPAS studies, the dissipation decreases with ESMF, GMLS-CAAS, WLS-ENOR, and TempestRemap schemes in that order.

Figure 11 L_min,2 metric for the TOPO field on RLL-CS meshes with different resolutions.

One of the key outcomes of this analysis indicates that using nonlinear filtering algorithms like CAAS can provide the benefits of property preservation to achieve global conservation and monotonicity constraints. However, these advantages can be offset by the higher local dissipation effects, especially on disparate mesh resolutions of Ω_s and Ω_t. So depending on the problem use case for which linear maps are being used to transfer solution fields between meshes with varying resolutions, an adaptive approach may be taken to choose a nonlinear filter when monotonicity constraints are important and to apply it with care, especially for high-order maps on which feature preservation may suffer.

5.5 Analysis on RRMs

We restrict the analysis in this study to all the conservative variations of the schemes chosen. As a result, the global conservation metric is truly satisfied for all cases tested with RRMs. The following subsections provide a detailed comparison of the error resolutions as well as global and local bound preservation metrics.

5.5.1 Error metric comparison

In order to compare the performance of the different remapping schemes on more realistic and complex unstructured RRMs, the global error norms ${∥E∥}_{L_{\mathrm{1}}}$, ${∥E∥}_{L_{\mathrm{2}}}$, ${∥E∥}_{L_{\mathrm{\infty }}}$, ${∥E∥}_{H_{\mathrm{1}}}$, and ${\left|E\right|}_{H_{\mathrm{1}}}$ are computed and tabulated in Table 7. These results were computed on the CS(r2)–MPAS(r2) RRM combination to observe the accurate resolution of different remapped fields.

We note that the global errors with respect to all error norms are considerably smaller in WLS-ENOR and TempestRemap for the smooth (analytical) field variables sampled on the RRM input meshes. Interestingly, the low-order ESMF implementation shows performance comparable to the high-order GMLS-CAAS(p=4) scheme in ${∥E∥}_{L_{\mathrm{2}}}$, even though both these methods are more than 5 orders of magnitude worse than WLS-ENOR in terms of absolute error value. The gradient error metrics (${∥E∥}_{H_{\mathrm{1}}}$ and ${\left|E\right|}_{H_{\mathrm{1}}}$) also show the high dissipative errors in low-order ESMF methods for these smooth field data.

However, it is imperative to note that when transferring field data with significant C¹ discontinuities (TOPO, CFR, TPW) that are typical in real-world climate simulations, using high-order remapping schemes such as WLS-ENOR, GMLS meshless methods and even TempestRemap(p=3) only provide nominal improvements over linear maps computed with low-order ESMF reconstructions with the conserve2nd option. These results demonstrate that the state-of-the-art ESMF conservative remapping schemes can be complemented by other algorithms that have relatively better accuracy profiles and lower feature dissipation for many coupled climate modeling problems.

Table 7 Comparison of global error norm metrics for different remap schemes on the finest (r2) CS-MPAS RRM combination, ordered from best to worst in terms of ${∥E∥}_{L_{\mathrm{2}}}$ error for all field variables.

In order to make a fair selection of ESMF conservative algorithm, the comparison between the conserve and conserve2nd implementations was repeated on RRMs. The ${∥E∥}_{L_{\mathrm{2}}}$ error profiles for various equal-resolution source and target meshes are shown in Fig. 12 for the TPW variable. The results suggest that the behavior of the conserve2nd option remains marginally better than the conserve option in these experiments, even though the magnitude of accuracy gain with conserve2nd is not as attractive in comparison to the case of regularly refined meshes analyzed in Fig. 6. Hence, as an outcome of this experiment, the recommendation is to use conserve2nd in all existing production climate simulations when high-order maps are unavailable, as it has demonstrated more favorable properties than the conserve option in ESMF.

Figure 12 ESMF comparison of conserve and conserve2nd options for CS-MPAS RRMs for the TOPO field using log(${∥E∥}_{L_{\mathrm{2}}}$) global error metrics as a function of N_r.

5.5.2 Monotonicity metrics

Similar to the uniformly refined cases, the monotonicity metric $\left|G_{\max }\right|$ was analyzed on the CS(r2)–MPAS(r2) regionally refined meshes for the CFR field, as plotted in Fig. 13. The behavior of the $\left|G_{\max }\right|$ metric shows similar trends in the RRM cases compared to the uniform refinement cases. The high-order TempestRemap method shows the worst behavior, while low-order ESMF maps exhibit good dampening after the first remap step. Note that the hidden traces of $\left|G_{\max }\right|$ for GMLS-CAAS and WLS-ENOR in Fig. 13 are at the expected levels of zero (monotone reconstructions). The observation shown here confirms that overlay-mesh-based schemes like ESMF and TempestRemap require additional post-processing monotone filters such as CAAS during the online solution projection step in a simulation in order to ensure strict global bounds preservation and to avoid Gibbs oscillations for strong gradient fields.

Figure 13 $\left|G_{\max }\right|$ metric for the CFR field on CS-MPAS RRM with (a) coarse to fine and (b) fine to fine resolutions.

5.5.3 Locality metrics

The observed trends in the locality metrics computed on RRMs are similar to those exhibited in the uniformly refined mesh experiments presented earlier. Hence, for brevity, only the $L_{max,\mathrm{2}}$ and L_min,2 metrics are presented for the TOPO field on the finest resolution (r2) of RRMs in Fig. 14.

It is imperative to recognize that the high-order remapping schemes from TempestRemap and WLS-ENOR continue to generate minimally diffusive projections on target RRMs. This behavior is quite attractive as local features in the fields, even in the presence of strong gradients, are resolved accurately with very low dissipation away from the sampled reference field data. The locally adaptive, discontinuity tracking WLS-ENOR algorithm shows strongly stable behavior for all test cases.

However, the meshless GMLS-CAAS scheme and the low-order ESMF scheme exhibit severe departures away from these local bounds that are consistent with the relatively larger ${∥E∥}_{H_{\mathrm{1}}}$ and ${\left|E\right|}_{H_{\mathrm{1}}}$ observed in Table 7. Depending on the solution fields being remapped between model components, such strong dissipation in sharp features could lead to local numerical artifacts in the coupled climate system.

One key observation from the global monotonicity metrics shown in Fig. 13 is that TempestRemap fails to maintain strict global bounds, while low-order ESMF and the meshless GMLS-CAAS schemes recover them consistently. However, the local bound preservation metrics in Fig. 14 demonstrate that when using the conservative, high-order TempestRemap algorithm, element-wise dissipation in the field data is strictly bounded and relatively much smaller in comparison to other remapping algorithms.

Figure 14 $L_{max,\mathrm{2}}$ (a) and L_min,2 (b) metric evolution for the TOPO field on the CS(r2)–MPAS(r2) RRMs.

6 Future research directions

This work provides a foundation for the systematic computation of key metrics of interest in regridding problems and implements the infrastructure to sample analytical and real-world fields as well as to measure the metrics through a flexible Python code. In the current paper, we have presented detailed experiments to compare remapping schemes for FV–FV discretization settings on three specific mesh families (CS, MPAS, RLL) of varying resolutions, both uniform and regionally refined. There are several research directions that will be explored further in order to strengthen and improve the intercomparison protocol and apply it to more realistic climate science problems.

Four key extensions of interest in this context are

1. complex meshes, including successive r-adapted meshes, and realistic samples with topological holes;

2. complex discretizations, other than just FV for source and target components;

3. complex field descriptions, including vector data and preservation of nonlinear correlations in multiple fields; and

4. computational efficacy, comparing numerical accuracy and algorithmic performance on next-generation architectures.

Further details regarding these research directions are provided in the following subsections.

6.1 Complex meshes: more than just global meshes

Uniformly refined meshes provide a simple infrastructure to test the consistency of remapping schemes, and additionally, regionally refined global grids help identify any spatial-scale dependencies in the algorithms. However, it is important to note that the meshes tested do not include any topological holes that are typical in production climate simulations. A direct extension of the current study would be to perform FV–FV remapping studies on realistic meshes, especially with the MPAS ocean polygonal meshes (Petersen, 2018) that do not cover the entire sphere.

Several potential issues may arise when remapping fields on these complex topological meshes: for example, dealing with narrow isthmus regions like Panama, which is also at the boundary between ocean basins, or along the coast of Peru where sea level with high-altitude points from the Andes can influence and contaminate remaps to produce biases that appear along these coastal regions. Such regions can yield incorrect field remaps when using simple nearest neighbor maps, high-order reconstructions, or even when computing the overlay mesh, as the numerical tolerances for querying neighborhood data become an important influencing factor in the overall accuracy.

The metrics definitions for standard error norms, global maxima, and global minima introduced in Sect. 3 are valid for general meshes. But, in order for the global conservation metric L_g shown in Eq. (11) to be valid, the target domain has to fully envelop the source domain. Additionally, the definition and implementation of local dissipation metrics use a small neighborhood (one-ring patch) around the target cell, which would also impose a similar requirement. This implies that for analysis with realistic meshes, the metrics for evaluating the quality of solution field data transferred between the atmosphere and ocean components will strictly be measured only on the global atmosphere meshes to satisfy the above constraints. Note that the input specification for such cases would require at least the area fractions and appropriate masks to be provided in order to enable efficient remap computations using the meshless and hybrid schemes like GMLS and WLS-ENOR.

Furthermore, another test study of interest would be to develop a sequence of meshes that slowly deform from the original mesh into some intermediary that resolves moving solution fields like tracer transport, eventually returning back to the original mesh. Such a test is referred to as “cyclic-remap” test suites in the ALE community (Bochev and Shashkov, 2005). Performing successive remap over the sequence of such meshes, one can directly compare the final resulting field after multiple remaps directly with the initial field without a need for a reference analytical solution field. Devising such a study requires mesh-based adaptivity and smoothing algorithms to be used effectively in addition to mesh optimization strategies to avoid element inversions (Brewer et al., 2003). On the other hand, any errors introduced due to SPH sampling functionals and discontinuities introduced due to clipping field bounds can be completely avoided with such a setup.

6.2 Complex discretizations: more than just finite volume

The current study focused primarily on element-averaged field data typical of FV discretization of models. Although many atmosphere and ocean models define the scalar fields as element-averaged piecewise constant data, it is imperative to extend the above analysis framework to spectral element (SE) data defined on CS meshes so that it can be remapped onto either FV data on MPAS or RLL meshes or SE data on CS meshes with different resolutions. These extensions would allow us to verify the flexibility of the high-order, conservative remapping schemes tested in the current study under discretization specifications commonly used in climate components like HOMME (Taylor et al., 2007). We also make a note that such extensions to more complex discretizations of source and target data should not require any modifications in the metrics definitions already presented in this paper. Further extensions to meshes that contain topological holes will still require Ω_s⊂Ω_t.

6.3 Complex fields: more than just scalars

During the remapping of climate data in production runs and scientific analysis, it is important to preserve not only scalar fields but also vector fields and derived properties that provide better insights from simulations. For example, to understand the atmospheric flows, analysis of global wind patterns is often performed in weather prediction, which requires the evaluation of derivatives or integrals of the vector wind velocity fields (Nair and Jablonowski, 2008). Some critical considerations in such scenarios require preservation of divergence-free velocity fields, conservation of vorticity, and preservation of wind-stress curl fields during the remap phase. Another example is, during the computation of tracer advection, certain tracers that contribute to an equation of state must be remapped consistently with the density field in order to preserve derived quantities.

Typically, such vector fields are treated as collections of unrelated scalar fields that are remapped independently. However, such approaches are deficient in preserving divergence-free conditions and can be inconsistent, since the propagation of remapping error in the components is not correlated. Care should also be taken to ensure that the regridding of vector components uses a proper 3D Cartesian coordinate system instead of spherical mesh projections, which will not yield consistent vector fields on Ω_t.

Remapping algorithms based on mimetic schemes (Pletzer and Hayek, 2019) that provide exact conservation for both scalar and vector fields are promising in this direction. To our knowledge, existing remapping algorithms based on L² minimization and high-order interpolation-based algorithms need further research for tackling vector field data in problems of interest. Furthermore, field data with nonlinear correlations such as those analyzed by Lauritzen and Thuburn (2012), adapted for remapping scenarios, can also be valuable in determining whether the solution transfer algorithms can remain conservative and preserve correlation properties. Nonlinear remapping schemes (Carey et al., 2001; Bochev and Shashkov, 2005; Bochev et al., 2011) may also prove to be viable options for these cases, although computational cost can be relatively much higher than using linear maps on decoupled scalar components. Additionally, some remapping metrics definitions introduced in Sect. 3 will have to be extended for vector field data.

6.4 Computational efficacy: comparing accuracy and efficiency

The additional characteristic dimension that is essential to recognize in intercomparison studies is that the overall cost to obtain the remaps over the simulation cycle includes both a one-time setup cost and a constantly growing computational load at every coupling step when field data need to be transferred between components. Hence, in order to better describe the computational cost of the remapping algorithms, we could split the effort into an offline cost and an online cost. Note that, often, what is sacrificed in terms of computational performance is gained in the quality of the remapped solution through higher accuracy, global bounds preservation, and strong local feature resolution with minimal dissipation. Since the numerical efficiency and time to compute the solution are competing factors, the usage of advanced remapping algorithms for realistic cases will require a more detailed analysis of the computational complexity at scale. An exception to the performance model derived here is that for dynamic or moving meshes such as those encountered in sea ice and ocean coupling in ESMs, the offline cost is nearly zero and the entire computational cost is online, since the maps to compute field projections have to be recomputed with changes in Ω_s and Ω_t.

The traditional mesh-based conservative remapping algorithms used in ESMF and TempestRemap require computation of an intersection or supermesh, Ω_s∩t, which in general is computationally expensive. While linear complexity strategies like advancing front methods do exist in libraries like PANG (Gander and Japhet, 2013) and MOAB (Mahadevan et al., 2020) to compute mesh intersections, they can suffer from robustness issues when Ω_s has topological holes. An alternative approach is to utilize a K–d tree data structure to accelerate queries on the meshes in order to exactly locate the element containing a particular point, which is a fundamental operation in intersection mesh computation (Jansen et al., 1992). The leading computational cost for this operation is 𝒪(Nlog (N)), where N is the total number of elements in the mesh. The linear map computation itself follows typical FV or FE operator assembly workflows and remains bounded. The online cost of such a linear map application is essentially then provided by the sparse matrix–vector (SpMV) operation. The complexity for SpMV is, however, dependent on the degree p used for field reconstruction as it dictates the bandwidth of coupling between degrees of freedom in Ω_s and Ω_t. The illustration in Fig. 15 for ESMF (conserve) and TR (p=3) remapping algorithms shows the number of CS elements required to reconstruct the solution on a single target MPAS element when computing projections between the CS(r3)–MPAS(r3) uniform refinement case. While the low-order ESMF algorithm may require only 2 FLOPS per reconstruction on a target element, TempestRemap (p=3) would require a minimum of 50 FLOPS to achieve better accuracy. To overcome some of these computational issues for large meshes, the use of MPI parallelism has been exploited in ESMF and has been proven to be scalable with MOAB–TempestRemap (Mahadevan et al., 2020) libraries.

Figure 15 Comparison of source field reconstructions showing DoF coupling in linear maps when utilizing different remapping algorithms and degrees of expansion.

Computation for the GMLS-CAAS mesh-free approach on a manifold is dominated by the internal QR with pivoting factorization, requiring $\mathcal{O}\left(\frac{(p^{\mathrm{2}}{)}^{\mathrm{3}}}{\mathrm{6}}\cdot N\right)$ FLOPS for the offline stencil calculation, whereby p is the polynomial degree of the basis used. This stencil calculation can be stored and applied as an SpMV operation, similar to TempestRemap. The CAAS nonlinear filter is an online cost as it is solution-dependent, so at every instance that it is applied it carries a $\mathcal{O}\left(\frac{p^{\mathrm{2}}}{\mathrm{2}}\right)$ bounds calculation cost for each of N degrees of freedom (DoF).

In terms of the computational cost, similar to TempestRemap and GMLS, the transfer for smooth function and the detection of discontinuities in WLS-ENOR can take advantage of some pre-processing steps to build matrix-based operators. Hence, they primarily involve SpMV as the primary online cost as well. Additionally, the resolution of discontinuities requires constructing and solving the generalized Vandermonde systems for each target cell near discontinuities as described in Sect. 2.2, which can be relatively expensive. Hence, the overall cost for resolving discontinuities depends on the percentage of the cells in discontinuity regions, which is determined only at runtime to account for field evolution. In general, on coarse meshes, this ratio may be relatively higher. As the underlying mesh is refined, the solution field profile has a finer approximation that results in a lower proportion of discontinuous cells. So the overall cost for the most expensive part in WLS-ENOR is expected to have a sublinear time complexity under mesh refinement.

While the theoretical complexity requirements for the remapping algorithms under consideration can be considerably different, it is imperative to measure the performance of these implementations on various architectures on standard test problems to gauge overall efficacy (accuracy vs. total computational time for N_r steps). This requires configuring, building, and installing the ESMF, TempestRemap, Compadre (GMLS-CAAS), and WLS-ENOR libraries on the same machine to compare the performance profile for test problems at scale. Given the complexity and magnitude of such a task, we express this as another future avenue for research experiments that can add value to the broader climate science community using remapping algorithms (Valcke et al., 2022).

7 Conclusions

Remapping algorithms are critically important in climate science applications to maintain key numerical properties during the transformation and transfer of coupled field data between component models. Inaccurate, nonconservative, or highly dissipative remapping operators can introduce numerical artifacts in the coupled fields, which can destroy the high-order accuracy of component solvers and propagate errors into the global nonlinear system representing the climate system. Hence, understanding the behavior of remapping algorithms for Earth system modeling requires standardized numerical definitions that provide better insight into the accuracy of transferred fields, along with the ability to preserve global and local solution bounds to avoid numerically induced instabilities in the system.

In this paper, with these motivations in mind, four different remapping schemes were selected based on the current availability of the software implementations, the maturity of the underlying numerics, and potential computational efficiency that could be obtained for real-world scenarios in comparison to existing state-of-the-art remapping algorithms used in production climate simulations. The mesh-based implementations in the libraries ESMF (Hill et al., 2004) and TempestRemap (Ullrich and Taylor, 2015; Ullrich et al., 2016) provide support for the projection of field data between unstructured grids by using intersection meshes to compute conservative weights through global L²-minimization approaches. These remapper implementations are used routinely in several Earth system solvers such as CESM (Hurrell et al., 2013) and E3SM (E3SM Project, 2018) for remapping scalar and flux field variables between component models. In contrast to the intersection-mesh-based remapping approaches, a hybrid mesh-based approach without an explicit need for overlay meshes using the WLS-ENOR algorithm (Li et al., 2020) was included in this study. Finally, a fully meshless remapping algorithm was selected, which is based on the generalized moving least squares (GMLS) scheme (Trask and Kuberry, 2020), with optional post-processing filters to ensure monotonicity in the projected solutions using the Clip-And-Assured-Sum (CAAS) (Bradley et al., 2019) technique. These remapping algorithms cover a large span of low- and high-order conservative solution transfer implementations that can directly impact the overall stability and accuracy of predictive solvers and analysis suites for weather and climate modeling.

Comparing these four distinct remapping algorithms requires a uniform test infrastructure, which provides the framework to create new verification studies and analyze the metrics obtained from the remapping algorithms. To enable such unbiased comparative studies, we introduced several remapping metrics that represent the key properties of remapping algorithms. These include global error measures under various norms like ${∥E∥}_{L_{\mathrm{1}}}$, ${∥E∥}_{L_{\mathrm{2}}}$, and ${∥E∥}_{L_{\mathrm{\infty }}}$ and gradient error measures given by ${∥E∥}_{H_{\mathrm{1}}}$ and ${\left|E\right|}_{H_{\mathrm{1}}}$. These error norms provide the necessary verification of the theory and implementation of the remapping algorithms by measuring the theoretical order of accuracy applied to solution fields with sufficient smoothness. Next, global conservation errors are measured by comparing the integral of the sampled reference fields on Ω_t against the remapped solution on Ω_t for multiple projection iterations. Finally, metrics for the global and local departure away from maxima or minima in the solution fields can be evaluated in various norms to provide insight on monotonicity preservation and feature dissipation due to the remapping algorithms. These standardized metrics provide the blueprint to build the remapping intercomparison suite, which was then used to understand the numerical properties of all algorithms under consideration.

Furthermore, a flexible workflow built on several Python-based drivers to generate the unstructured meshes of different element topologies and resolutions, including regionally refined meshes, and to accurately sample five element-averaged fields using SPH expansions was provided. With these input meshes, the four remapping algorithms were applied for both the smooth analytical fields and representative real fields in an iterative fashion to compute cyclic projections (${\mathrm{\Omega }}_{\mathrm{s}}\to {\mathrm{\Omega }}_{\mathrm{t}}\to {\mathrm{\Omega }}_{\mathrm{s}}$) for FV–FV field transfers between component models in a climate system.

(Guerra et al., 2021)(Guerra et al., 2021)Guerra and Mahadevan (2021)(Wimmers and Velden, 2011)(Platnick et al., 2020)(NOAA National Geophysical Data Center, 2009; Amante and Eakins, 2009)(Guerra and Mahadevan, 2021)Mahadevan et al. (2021)(Mahadevan et al., 2021)

Table 8 Code and data availability.

The results compiled from various test problems demonstrate that the conservative remapping implementations in ESMF are all first-order accurate for smooth problems, even though the conserve2nd option produces consistently better accuracy than the conserve option. Hence, when possible, the ESMF conserve2nd option should be used to obtain better remaps for climate simulations. However, these low-order ESMF maps are highly dissipative in general. In contrast, TempestRemap produces remaps that are globally conservative with high-order accurate convergence to 𝒪(h^p+1) for smooth solution fields using second-degree (p=2) polynomial reconstructions and bounded dissipation on field projections containing sharp features. Here, h represents the characteristic spatial length of the mesh wherein errors are measured. However, neither of these low- and high-order L²-minimization approaches can guarantee monotone remaps without the use of post-processing filters such as slope limiters or CAAS to enforce global and local solution bounds. The use of these filters, however, inevitably introduces additional dissipation that can become significant with repeated applications of the remap operator.

On the other hand, the hybrid mesh-based and meshless schemes achieve very high order consistently (up to 𝒪(h⁵)) while retaining global conservation, which has been one of the key advantages of traditional overlap-mesh-based schemes implemented in ESMF and TempestRemap. Additionally, the ability to capture smoothly varying fields very accurately without any degradation even for high-polynomial-degree reconstructions makes these schemes attractive and competitive compared to traditional mesh-based schemes utilizing L²-minimization methods for ESMs that require computation of overlap meshes, which can incur a significant computational cost at high spatial resolutions (Mahadevan et al., 2020). These methods can be especially valuable when highly accurate scalar fields need to be sampled on a refined RLL grid for further analysis or in situ visualization to track multidecadal climate evolution.

It is essential to note that while the fully mesh-based remapping algorithms are, in general, insensitive to mesh resolutions or the topology, both the hybrid and meshless schemes are susceptible to larger mesh-dependent dissipation. This is especially evident when the source and target mesh resolutions differ drastically. The use of the CAAS algorithm, when combined with even nonconservative schemes such as GMLS, can provide global conservation and bounds preservation at the cost of added dispersion to control numerical oscillations. In contrast, the built-in discontinuity indicators used by the WLS-ENOR algorithm demonstrate good feature-resolving properties for real fields in all experiments conducted.

These experiments conducted on both uniform-resolution meshes and regionally refined meshes provide valuable insight into the properties of remapping algorithms and their numerical behavior. However, practical use of these algorithms in real-world scenarios requires deeper investigation using more topologically diverse, complex meshes and discretization specifications that are more representative of components used in E3SM and CESM.

Finally, we want to emphasize that the MIRA infrastructure presented in this paper is freely available as an open-source package (Guerra et al., 2021) to compare new and existing remapping algorithms under the same overall test constraints. Such intercomparison studies are important to evaluate the cost of remapping algorithms under stability and accuracy constraints, which remain crucial to better understanding the propagation of errors in coupled climate and weather systems.