Metrics for Intercomparison of Remapping Algorithms (MIRA) applied to Earth System Models

Strongly coupled nonlinear phenomena such as those described by Earth System Models (ESM) are composed of multiple component models with independent mesh topologies and scalable numerical solvers. A common operation in ESM is to remap or interpolate results from one component’s computational mesh to another, 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 ESM. However, a unified approach to compare the properties of these different schemes has not been attempted previously. We present 5 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 are 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), 10 TempestRemap (Ullrich and Taylor, 2015), Generalized Moving-Least-Squares (GMLS) (Trask and Kuberry, 2020) with postprocessing 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 the 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 15 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 highorder 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 20 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 1 https://doi.org/10.5194/gmd-2021-323 Preprint. Discussion started: 28 October 2021 c © Author(s) 2021. CC BY 4.0 License.

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 are 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), 10 TempestRemap (Ullrich and Taylor, 2015), Generalized Moving-Least-Squares (GMLS) (Trask and Kuberry, 2020) with postprocessing 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 the 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 15 remapping performance when dealing with non-smooth data. There is a failure to maintain monotonicity in the traditional L 2 -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 highorder 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 20 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 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 2 -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 5 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 10 mesh, field, and scale resolution-independent metrics definitions. This intercomparison paper is organized as follows: -A detailed literature survey on the current state-of-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 Section 2.
-Definitions for remapping metrics that evaluate field accuracy, global conservation, strict global bounds control, and 15 feature dispersion are presented in Section 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 Section 4. Some implementation details on the open-source Python-based infrastructure are also provided. 20 -Consolidated results from the intercomparison study applied to four competing remap algorithms on representative problems are shown in Section 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 Section 6, followed by a summary of the intercomparison study in Section 7.
2 Remapping Algorithms 25 In general, for coupled simulations, we need to transfer a field ψ s ∈ R ns defined on the source grid, Ω s , to a quantity ψ t ∈ R nt on the target grid, Ω t . A remap is hence defined as an operator R : R ns → R nt that transfers ψ s to a trial field ψ t ∈ R nt . As noted above, we generally wish to preserve a property defined as an operator P : R nt → 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 expansion using nodal basis functions.

Related Work
One of the simplest data-transfer methods is to use piecewise interpolation functionals. This approach is particularly convenient 5 if ψ s defined on Ω s has an associated function space, which can be directly utilized for evaluating the interpolator. Such consistent interpolation techniques utilizing the underlying basis functions for field descriptions give rise to standard secondorder 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. 10 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 15 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 its 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 reconstructing locally the manifold geometry from a point set representation, and then generating a compact stencil in the local coordinate 20 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 C 0 25 normal field, instead of global stencils and C ∞ 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); Jiao and Heath (2004a). Note that the high-order quasi-interpolants (with degree of 30 basis expansion p > 1) tend to be significantly less dissipative than low-order methods (p ∈ [0, 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, 5 https://doi.org/10.5194/gmd-2021-323 Preprint. Discussion started: 28 October 2021 c Author(s) 2021. CC BY 4.0 License. 1994;Grandy, 1999;Gander and Japhet, 2013;Jones, 1999;Ullrich and Taylor, 2015) that rely on common-refinement-based L 2 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 5 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 (aka C 0 discontinuities), which 10 tend to lead to O(1) oscillations (or "ringing") that do not vanish under mesh refinement, analogous to the Gibbs phenomena (Gottlieb and Shu, 1997;Fornberg and Flyer, 2007). Additionally, any discontinuities in the derivatives (C 1 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 non-physical variations that can influence the numerical stability of the coupled global nonlinear multiphysics system. 15 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, often the properties of the linear maps can be enhanced by using a procedure that is nonlinear (depending on the 20 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 25 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 on 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 30 past, a discontinuity detecting, a posteriori stabilization procedure has 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 5 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 is 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 10 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.

15
A common theme across all of the remapping methods described in this paper is that they utilize some variants of the least squares approximations (aka quasi-interpolation) in their computational kernels. To illustrate the idea, let us consider a function f (u) : Ω → R at a given point u 0 = [0, 0] T , such as a quadrature point in a cell on the target mesh. Let us first assume a domain Ω ⊂ R 2 for simplicity, where u ≡ [u, v]. Let f be C p−1 continuous for some degree p ≥ 0, then f (u) can be approximated to p + 1 st order accuracy about u 0 , using a degree-p two-dimensional Taylor polynomial as where c jk = 1 , and h is a measure of the radius of the local neighborhood. In cell-to-cell transfer, typically some integrals of f are known on the source mesh. Given m cells on the source mesh in a neighborhood of a point u 0 on the target mesh, and 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 + 1)(p + 2)/2 unknown coefficients c jk , for 1 ≤ i ≤ m. Note that one can also use a bi-degree-p Taylor polynomial by letting 0 ≤ j ≤ p and 0 ≤ k ≤ p in (1), which would lead to n = (p + 1)(p + 1) unknowns. The resulting approximate Taylor polynomial can then be used, for example, to evaluate (or quasi-interpolate) f at u 0 to p + 1 st 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 j v k with x and x j y k z , respectively. The m equations in (2) can be written in matrix form as Ax ≈ b, where A ∈ R m×n is known as a generalized Vandermonde matrix, x ∈ R n is composed of c jk in (1), and b ∈ R 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 r = b − Ax, i.e., where W is a diagonal matrix containing the weight for each row in A. This formulation is known as weighted least squares or WLS in 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 10 the weighted least-squares problem in Equation (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).

Algorithms for Earth System Models
Among the standard conservative remapping and high-order reconstruction strategies introduced in Section 2.1 for climate 15 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 ESM's, 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 2 projection support Both ESMF Regrid and TempestRemap provide the remapping capability for scalar fields defined on Ω s based on a 25 weighted residual formulation given in Equation (3). With a source field ψ s defined on Ω s , the formulation computes ψ t on Ω t by solving the following problem: where the φ i are suitable weight functions. In a common-refinement-based transfer (Jiao and Heath, 2004a), the φ i are 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 Equation (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 does provide both first-and second-order conservative remapping 5 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 method (GMLS), with global 10 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 15 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 Equation (3).
Other opportunities for comparison in this category include using reconstructed climate data with the conservative bilin-20 ear algorithm in ESMF or bicubic interpolations available in YAC (Hanke et al., 2016).
These algorithms and their implementations span a range of remapping techniques, including those used currently 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.

25
The Earth System Modeling Framework (ESMF * ) (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 modeling. 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 non-conservative nearest neighbor interpolants. ESMF can produce maps that are "offline", in the sense that 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 5 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 10 of the algorithms in ESMF, we refer readers to Collins et al. (2005), Jones (1999), and Kritsikis et al. (2017).

TempestRemap
TempestRemap † is a software library to generate conservative, consistent, and monotone remapping weights of arbitrarily high order accuracy 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 15 volume and/or spectral element (both continuous and discontinuous spaces) descriptions with conservative prescriptions. For purposes of the current manuscript, 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 20 are computed according to a weighted pseudoinverse 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 Ω s∩t = Ω s ∩ Ω 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 Equation (2) can be very high, leading to numerical roundoff errors and degradation in the overall accuracy of the remap 25 operator to first-order accuracy in the neighborhood.
A common approach for the integration operator has been to transform integrals over mesh faces into line integrals over the boundary using the divergence theorem (Dukowicz and Kodis, 1987). 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 30 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. † TempestRemap: https://github.com/ClimateGlobalChange/tempestremap To obtain conservation, the coefficients of the matrix are projected linearly into the space of maps that are both conservative and consistent.

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, 5 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' 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 10 approximation and wellposedness. 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 (H(div)), and Nedelec basis (H(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 15 in these various non-standard 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 When a GMLS problem is solved, a set of coefficients corresponding to elements of a basis (traditionally polynomial) are computed. The target operation is then applied to each member of the basis used for approximation, after which a dot product 20 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 25 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 3 in R 2 , as compared to a traditional basis in R 3 , 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-30 dependent and generally vanishes with refinement. However, such weak conservation notions may not be deemed satisfactory for climate modeling applications where 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, global bounds, along with an attempt to improve local property preservation.
Similar to the overlay-mesh-based, ESMF, and TempestRemap offline remappers described in Section 2.3.1 and Section 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 between meshes, and in the process, it enhanced the treatment of discontinuities to resolve not only the severe oscillations at C 0 discontinuities (i.e., jumps in function values), but also the accumulated effect of mild oscillations due to C 1 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 30 if available). More specifically, WLS-ENOR utilizes adaptive quadrature rules with p-refinement in smooth regions and hrefinement near discontinuities in its numerical integration to achieve high-order accuracy and (near-local) conservation. More specifically, WLS-ENOR is composed of three major components. The first component is a WLS-based algorithm for smooth ‡ Compadre Toolkit: https://github.com/SNLComputation/compadre functions, as we described in Section 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 cancellations 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 inter-comparison 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 (aka Heaviside functions) over 5 the cells on the source mesh as the test functions in (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 2-D object embedded in R 3 . One could construct a WLS in R 3 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 0 ∈ S. Specifically, let m 0 denote an approximate normal at x 0 , which can be the exact normal to S or a first-order estimation.

10
Let t 1 and t 2 form an orthonormal basis of an approximate tangent plane orthogonal to m 0 . The local uv coordinate frame is then centered at x 0 with axes t 1 and t 2 . We can then transform the sample points about x 0 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.

15
The second component in WLS-ENOR is the detection and resolution of discontinuities. In particular, WLS-ENOR can detect C 0 and C 1 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 C 0 and C 1 discontinuities as described in (Li et al., 2020). The original detector in (Li et al., 2020), however, was for node-tonode transfer. For cell-to-cell transfer in this work, we simply apply the detector on the dual mesh and treat cell-averaged 20 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 (cf. Section 2.1). In terms of implementation, the discontinuity detector on the source mesh is implemented as a 25 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, 30 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 sub-divided 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 5 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, where 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 10 the future for both node-to-node and cell-to-cell field transfers.

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 15 metrics can be grouped. 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 section 4) for the realistic fields. The regridding operator from 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 25 these operators take the form where J j denotes the area or weight of the j th degree of freedom (i.e. the volume of a finite volume).

Grid Sensitivity
A crucial factor for the success and broader usability of a general remapping algorithm in ESMs is the ability to produce meshindependent 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

Standard Accuracy Measures
Accuracy in the remapped solution will be assessed with standard error metrics defined as follows.
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 Equation (5) on Ω t . In some general sense, the error measure E L1 identifies errors in large-scale features, E L2 identifies errors in small-scale features, and E L∞ identifies the largest pointwise error. Related to accuracy, consistency is assessed by applying these metrics to uniformly smooth fields with no C 0 or C 1 discontinuities, and verifying the asymptotic theoretical rate of convergence under 10 uniform mesh refinement conditions.
Note that in order to eliminate potential aliasing errors, the normalization factors D t [ψ] used in the denominator for definitions of E L1 , E L2 , and E L∞ 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 .

15
Preservation of the solution gradients in addition to other critical properties, such as local conservation in the remapping pro- . Let ∇D s , ∇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: |E| H1 seminorm and the E H1 norm. where is the corresponding gradient of the error, and I t is the weighted integral using Equation (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.

Global conservation
Global conservation is trivially assessed by evaluating the change in the global integral of the scalar field value on the source 5 mesh and the projected field on the target mesh. We use the following metric to quantify global conservation.
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 Equation (11) would be for remapping heat and moisture fluxes from the land surface to the overlying atmosphere.

Global Extrema Preservation
Global extrema preservation can be assessed via the standard G min and G max error metrics (Ullrich and Taylor, 2015): The error measures G min and G max identify undershoots and overshoots, respectively, by taking on nonzero values when 15 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 presence of Gibbs phenomenon. Hence, the global extrema metric is particularly useful as it provides indications about the monotonicity-preserving properties of the remapping schemes.

Local Extrema Preservation
Local extrema preservation can be assessed using a localized difference, i.e. to what degree does the remapped grid cell value 20 fall within the range of surrounding grid cells sampled on the target grid. This consideration motivates us to define a localized difference in extrema: 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 : Note that the definition of the localized differences shown in Equation (14) and Equation (15) utilize a local neighborhood 5 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 a O(h) dependence on the mesh resolution and can be applied to C 0 or smoother fields, but not when C 0 discontinuities are present.

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 ) 10 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, uni-directional 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 Section 3. Here, N r indicates the number of iterative applications of the linear map to compute field transformation on Ω s → Ω t → Ω s , which will be referred to as remap iterations. The additional remapping 15 metrics computed from the repeated remap iterations are critical in determining the long-term temporal behavior of a fully coupled climate simulation, as it entails multiple transfers of field data between component models. Hence, determining the stability of the remapping operator and understanding dissipation behavior for repeated transfers provide valuable insight into the propagation of spatial coupling errors in the solver.

20
The workflow necessary to evaluate a given remapping method comprises of five consecutive steps described below. 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 5 fields by sampling from either an analytical function on the sphere or a set of prescribed Spherical-Harmonic (SPH) coefficients, which is described in Section 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 compute the metrics defined in Section 3 accurately.

10
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 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 output from each of the algorithms, 15 for the test problems devised, are the discrete solution vectors ψ i s ∈ Ω s , and ψ i t ∈ Ω t , where i ∈ [1, . . . N r ]. In the current study, unless otherwise specified, N r = 1000. The schematic shown in Figure 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 available publicly .

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 and the sampling methodology used in the Python implementation are provided below.

Idealized Analytical Fields
Idealized fields used in this study mirror the approach of Lauritzen et al. (2007) 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 Figure 3.
The first analytical field (ANALYTICALFUN1) is a combination of spherical harmonics functions with frequency wave 5 similar to order 3, given by where Y l m are the real spherical harmonic functions evaluated through the SHTOOLS package for degree m and polynomial order l.
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).

5
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 data set, 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 10 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 Figure 4 with a corresponding linear fit, we then produce controlled randomized realizations for each field on any unstructured mesh and 15 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 Figure 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 . The original data ranges from 0.25 to 0.1 degrees 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 20 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.  Global Topography (TOPO): Global topography data is taken from the ETOPO1 Global Relief Model || (Center, 2009;Amante and Eakins, 2009). This field has a global minimum/maximum for all realizations of -10994.0/8848.0 m, but is otherwise smoothly varying. A randomized reconstruction of TOPO is shown in Figure 5c.
The raw satellite data snapshot of TPW, CFR, and TOPO fields that are used in the current study have been made available separately  in order to make the workflow reproducible. , 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. 20 Detailed results from the intercomparison study and discussion on the implication of each metric to the remapping scheme are presented next.

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 25 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 O(h p+1 ) in E L1 , E L2 , and E L∞ global error norms, and O(h p ) in E H1 and |E| H1 global gradient error norms.
In these studies, we observed that convergence rates for both low and high-degree approximations up to p = 2 of the analyti- 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.

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 5 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 1. The ESMF first-order conservative scheme yields expected rates. However, the second-order scheme shows degraded convergence rates as confirmed by the global error and gradient norms. This convergence 10 result is unexpected and contrary to the analysis of the second-order, piece-wise linear finite-volume reconstruction procedure presented by (Kritsikis et al., 2017), which is implemented in ESMF. 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 Figure 6, where the E L2 and E H1 error profiles for both the conservative schemes are compared as a function 15 of remap iterations (N r ) for the CS-MPAS mesh combination. The computed error metrics from both the schemes indicate that, even though the convergence rates are similar, the 'conserve2nd' option in ESMF clearly 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.

TempestRemap
The conservative high-order linear maps computed by TempestRemap, as shown in Table 2, 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.

5
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∞ shows further reduction in convergence rates. The failure to achieve O(h 4 ) or higher-order accuracy is an artifact of the implementation in 10 TempestRemap, which can be improved with further analysis and is not a limitation of the underlying numerical method.  TempestRemap produces conservative solution projections between mesh combinations that can be 3rd-order with p = 2, for smooth fields. However, even if local element-wise conservation can be guaranteed in overlay-based high-order L 2 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 5 encountered, global bounds can be preserved in TempestRemap by employing nonlinear filtering algorithms such as CAAS during the online solution transfer in ESMs.

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 10 shown in Table 3 for various polynomial degrees.  The convergence rates for the nominal GMLS scheme and the GMLS-CAAS remapping method with a post-processing step in Table 3 show that, in general, high-order accuracy can be achieved. However, using the nominal GMLS scheme for climate modeling problems can result in non-conservative 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

WLS-ENOR
The convergence analysis for the WLS-ENOR scheme shows that even for high polynomial degrees, theoretically expected rates are observed for the smooth analytical field. The convergence rates for various degrees are tabulated in Table 4. 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 5 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.

Real Fields: Convergence Rate Comparisons
While high-order convergence rates are achievable for smooth field profiles of ANALYTICALFUN1 and ANALYTICALFUN2, 10 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 C 0 and C 1 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 5. The convergence rate comparisons using different remap schemes show that the computed rates for higher-order polynomial reconstructions fall severely short of theoretical rates 15 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 H1 and |E| H1 . The meshless schemes are competitive in terms of accuracy bounds and provide a viable alternative approach for usage in production climate simulations, where traditionally L 2 -minimization schemes and low-order bilinear maps have been routinely used. 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, typically p = 1, or p = 2 may yield a more numerically stable solution due to Gibbs phenomena. In such operational circumstances, high-order methods for TempestRemap and GMLS should be 5 augmented with nonlinear post-processing filters, such as CAAS.

Global Conservation
All mesh-based L 2 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 10 sum operations on the discrete linear map matrix representing the projection or mass-matrix operator.
The global field integral error indicating conservation deficit for TPW field variable is shown in Figure 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 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 Figure 7, have been verified for both low and higher-order maps, irrespective of the mesh topology, or the field variables being transferred across meshes.
Notably, the slow accumulation of roundoff error in the global field integral as a function of remap iterations has a detrimental effect to global conservation. This effect is seen in the field projections computed with TempestRemap through two sparse matrix-vector (SpMV) products in each remap iteration and is distinctly subdued in the low-order ESMF map projections, The WLS-ENOR scheme achieves excellent global conservation by applying adaptive quadrature rules to integrate the functions to (near) machine precision and redistributing the conservation errors locally near discontinuous regions. The unmodified GMLS scheme is non-conservative 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 a conservation constraint may not be strictly mandated, and this enables the use of even the non-conservative GMLS scheme, which has been demonstrated to achieve excellent accuracy.

5
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 Figure 8 and Figure 9 for L max of CFR variable on MPAS-RLL, and L min 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 10 violated.
The results clearly demonstrate that the mesh-based remapping schemes 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 recovered after several remaps (around 700 in Figure 8). However, the global minima increase 15 continuously within the bounds of the reference solution. Hence, the 'conserve2nd' option in ESMF damps the global maxima and amplifies the global minima.
In contrast, the high-order TempestRemap solutions show a drift away from the bounds in every case. The mildly damped increase of 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 20 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 and hybrid remap schemes, namely GMLS-CAAS and WLS-ENOR, maintain the global bounds for all test cases and fields tested in this study. Augmenting the non-monotone GMLS scheme with CAAS bounds preservation algorithm shows excellent, stable behavior for all cases. WLS-ENOR shows relatively 25 high compliance to 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 Figure 8(b).

30
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 com- These results showcase 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 5 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 its use of adaptive quadrature rules.

5
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 1-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. Lmax,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 higherorder 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 Figure 10 and Figure 11.  Figure 11. Lmin,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 5 care, especially for high-order maps where feature preservation may suffer.

Analysis on RRM Meshes
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 RRM meshes. The following subsections provide a detailed comparison of the error resolutions and global/local bound preservation metrics.

Error Metrics Comparison
In order to compare the performance of the different remapping schemes on more realistic and complex unstructured RRM meshes, the global error norms, E L1 , E L2 , E L∞ , E H1 and |E| H1 are computed and tabulated in Table 6. These results were computed on the finest CS-MPAS RRM mesh combination to look at the accurate resolution of different fields.
The global errors with respect to all error norms are considerably smaller in WLS-ENOR and TempestRemap for all the field 5 variables sampled on the RRM input meshes.
In contrast, the low-order ESMF implementation results in a larger error comparable to GMLS-CAAS with p = 4 for all the real fields, even though the relative accuracy on smooth analytical fields is much better with GMLS-CAAS. The use of ESMF with 'conserve2nd' option and the GMLS-CAAS schemes deliver larger E H1 and |E| H1 errors than the more accurate WLS-ENOR and TempestRemap remapped solutions. This summary shows that the state-of-art ESMF conservative 10 remapping schemes can be complemented by other algorithms that have relatively better accuracy profiles for many coupled climate modeling problems.
In order to make a fair selection of ESMF conservative algorithm, the comparison between the 'conserve' and 'conserve2nd' implementations was repeated on RRM meshes. The E L2 error profiles for various equal resolution source and target meshes are shown in Figure 12 for the TPW variable. The results suggest that the behavior of 'conserve2nd' option remains marginally 15 better than the 'conserve' option in these experiments, even though the magnitude of accuracy gain with the 'conserve2nd' is not as attractive in comparison to the case of regularly refined meshes analyzed in Figure 6. An outcome of this experiment is that the results clearly demonstrate the superiority of the 'conserve2nd' conservative remapping implementation option in ESMF in comparison to the 'conserve' option, and is hence recommended for use in all existing production climate simulations when high-order maps are unavailable.

Monotonicity Metrics
Similar to the uniformly refined cases, the monotonicity metric L max was analyzed on the finest resolutions of the CS-MPAS regionally refined meshes for the CFR field, as plotted in Figure 13.

Locality Metrics
The observed trends in the locality metrics computed on RRM meshes are similar to those exhibited in the uniformly refined mesh experiments presented earlier. Hence, for brevity, only the L max,2 , and L min,2 metrics are presented for the TOPO field on the finest resolution of RRM meshes in Figure 14.
It is imperative to recognize that the high order remapping schemes from TempestRemap and WLS-ENOR continue to 5 generate minimally diffusive projections on target RRM meshes. 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 departure away from these local bounds that are consistent with the relatively larger E H1 and |E| H1 observed in Table 6. Depending on the solution 10 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 Figure 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 Figure 14 demonstrate that when using the conservative, high-order TempestRemap algorithm, 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 10 3. Complex field descriptions, including vector data and preservation of nonlinear correlations in multiple fields 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.

Complex Meshes: More than just global meshes
Uniformly refined meshes provide a simple infrastructure to test the consistency of remapping schemes, and additionally, 15 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 5 these coastal regions. Such regions can yield incorrect field remaps when using simple nearest-neighbors 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 Section 3 are valid for general meshes. But, in order for the global conservation metric L g shown in Equation (11) to be valid, the target domain has 10 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 15 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 20 directly with the initial field without a need for a reference analytical solution field. Devising such a study requires meshbased 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

25
The current study focused primarily on element averaged field data that is 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 30 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 manuscript. Further extensions to meshes that contain topological holes will still require Ω s ⊂ Ω t .

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 5 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 10 components are not correlated.
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 2 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 15 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, albeit computational cost can be relatively much higher than using linear maps on decoupled scalar components. Additionally, some remapping metrics definitions introduced in Section 3 will have to be extended for vector field data. 20

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 25 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.
The traditional mesh-based conservative remapping algorithms used in ESMF and TempestRemap require computation of 30 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 Kd-tree datastructure 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 O (N log(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 does remain 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, 5 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 Figure 15 for  Computation for the GMLS-CAAS meshfree approach on a manifold is dominated by the internal QR with pivoting factorization, requiring O (p 2 ) 3 6 · N FLOPS for the offline stencil calculation, where p is the polynomial degree of the basis used. This stencil calculation can be stored and applied as a SpMV operation, similar to TempestRemap. The CAAS nonlinear filter https://doi.org/10.5194/gmd-2021-323 Preprint. Discussion started: 28 October 2021 c Author(s) 2021. CC BY 4.0 License.
is an online cost as it is solution dependent, and so at every instance that it is applied it carries a O( p 2 2 ) bounds calculation cost for each of N DoFs.
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 preprocessing 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 5 and solving the generalized Vandermonde systems for each target cell near discontinuities as described in Section 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 10 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 perfor- 15 mance 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.

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, non-conservative, or highly 20 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. 25 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-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 30 intersection meshes to compute conservative weights through global L 2 -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 remap-https://doi.org/10.5194/gmd-2021-323 Preprint. Discussion started: 28 October 2021 c Author(s) 2021. CC BY 4.0 License. ping 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 ClipAndAssuredSum (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 5 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 L1 , E L2 , E L∞ , and gradient error measures given by E H1 10 and |E| H1 . 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 15 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 20 analytical fields and representative real fields in an iterative fashion to compute cyclic projections (Ω s → Ω t → Ω s ) for FV-FV field transfers between component models in a climate system.
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 25 for climate simulations. In contrast, TempestRemap produces remaps that are globally conservative and high-order accurate up to O(h 3 ) convergence rates for smooth solution fields using p = 2. These low order ESMF maps are highly dissipative in general and can be rectified even with O(h 2 ) maps produced with TempestRemap, which show bounded dissipation even for sharp features. However, neither one of these low and high-order L 2 -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 30 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 O(h 5 )) 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, as compared to traditional mesh-based schemes utilizing L 2 -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 especially be 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. 5 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 liable 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 non-conservative 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 10 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. 15 Finally, we want to emphasize that the MIRA infrastructure presented in this paper is freely available as an open-source package  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. The Python intercomparison workflow infrastructure and the scripts to compute the metrics data for remapping schemes are available . This code was developed through funding from the CANGA project, and is publicly released under an open-source license, with copyright owned by UChicago Argonne, LLC. v1.0 was tagged in September 2021 and the archive is available for download from Zenodo. The original repository is hosted in GitHub. DOI: 10.5281/zenodo.5518037.

Raw satellite data snapshot
The raw satellite data snapshots used to generate the Spherical Harmonic decomposition for TPW, CFR, and TOPO fields are available at . The files are distributed under the Creative Commons Attribution 4.0 International License, while acknowledging the original sources for the data for TPW (Wimmers and Velden, 2011), CFR (Platnick et al., 2020), and TOPO (Center, 2009;Amante and Eakins, 2009). v1.0 was released on Aug 09, 2021 and available here. DOI: 10.5281/zenodo.5172792.

Input meshes and output metrics data
The pre-processed input meshes that were used in the study are available at . The input meshes and the consolidated metrics data for each of the remapping methods Author contributions. VM, JG, XJ and PK wrote the paper (with several helpful review comments from PB, PU and RJ). JG wrote significant portions of the remapping intercomparison code , with contributions from VM, and PK. VM generated the remap output data using ESMF and TempestRemap libraries. PK generated results using GMLS, and GMLS-CAAS. XJ and YL generated results using WLS-ENOR algorithms. VM consolidated the final metrics data from the remapping teams, performed analysis of the results, and generated the necessary data for the algorithmic intercomparison study. The broader project idea was conceived by PJ and PU.
Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-06CH11357, to generate the output for test problems using ESMF and TempestRemap remapping software libraries. We also gratefully acknowledge use of the Common Engineering Environment (CEE) Compute Servers at Sandia National Laboratories to generate the output for test problems using GMLS. Computational results for WLS-ENOR were obtained using the Seawulf cluster at the Institute for Advanced Computational Science of Stony Brook University, which was partially funded by the Empire State Development grant NYS #28451.