Metrics for Intercomparison of Remapping Algorithms (MIRA) protocol applied to Earth system models
 ^{1}Mathematics and Computational Science Division, Argonne National Laboratory, Lemont, IL 60439, USA
 ^{2}OU/CIMMS, NOAA National Severe Storms Laboratory, Norman, OK, USA
 ^{3}Department of Applied Mathematics & Statistics, Stony Brook University, Stony Brook, NY 11704, USA
 ^{4}Center for Computing Research, Sandia National Laboratories, P.O. Box 5800, Albuquerque, NM 87125, USA
 ^{5}Department of Land, Air and Water Resources, University of California, Davis, CA 95616, USA
 ^{6}Department of Mathematics, University of California, Davis, CA 95616, USA
 ^{7}Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
 ^{1}Mathematics and Computational Science Division, Argonne National Laboratory, Lemont, IL 60439, USA
 ^{2}OU/CIMMS, NOAA National Severe Storms Laboratory, Norman, OK, USA
 ^{3}Department of Applied Mathematics & Statistics, Stony Brook University, Stony Brook, NY 11704, USA
 ^{4}Center for Computing Research, Sandia National Laboratories, P.O. Box 5800, Albuquerque, NM 87125, USA
 ^{5}Department of Land, Air and Water Resources, University of California, Davis, CA 95616, USA
 ^{6}Department of Mathematics, University of California, Davis, CA 95616, USA
 ^{7}Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Correspondence: Vijay S. Mahadevan (mahadevan@anl.gov)
Hide author detailsCorrespondence: Vijay S. Mahadevan (mahadevan@anl.gov)
Strongly coupled nonlinear phenomena such as those described by Earth system models (ESMs) are composed of multiple component models with independent mesh topologies and scalable numerical solvers. A common operation in ESMs is to remap or interpolate component solution fields defined on their computational mesh to another mesh with a different combinatorial structure and decomposition, e.g., from the atmosphere to the ocean, during the temporal integration of the coupled system. Several remapping schemes are currently in use or available for ESMs. However, a unified approach to compare the properties of these different schemes has not been attempted previously. We present a rigorous methodology for the evaluation and intercomparison of remapping methods through an independently implemented suite of metrics that measure the ability of a method to adhere to constraints such as grid independence, monotonicity, global conservation, and local extrema or feature preservation. A comprehensive set of numerical evaluations is conducted based on a progression of scalar fields from idealized and smooth to more general climate data with strong discontinuities and strict bounds. We examine four remapping algorithms with distinct design approaches, namely ESMF Regrid (Hill et al., 2004), TempestRemap (Ullrich and Taylor, 2015), generalized moving least squares (GMLS) (Trask and Kuberry, 2020) with postprocessing filters, and WLSENOR (Li et al., 2020). By repeated iterative application of the highorder remapping methods to the test fields, we verify the accuracy of each scheme in terms of their observed convergence order for smooth data and determine the bounded error propagation using challenging, realistic field data on both uniform and regionally refined mesh cases. In addition to retaining highorder accuracy under idealized conditions, the methods also demonstrate robust remapping performance when dealing with nonsmooth 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 meshbased WLSENOR schemes. Local feature preservation analysis indicates that highorder methods perform better than loworder dissipative schemes for all test cases. The behavior of these remappers remains consistent when applied on regionally refined meshes, indicating meshinvariant implementations. The MIRA intercomparison protocol proposed in this paper and the detailed comparison of the four algorithms demonstrate that the new schemes, namely GMLS and WLSENOR, are competitive compared to standard conservative minimization methods requiring computation of mesh intersections. The work presented in this paper provides a foundation that can be extended to include complex field definitions, realistic mesh topologies, and spectral element discretizations, thereby allowing for a more complete analysis of productionready remapping packages.
The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”) and Sandia National Laboratories. Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DEAC0206CH11357, and Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DENA0003525.
The U.S. Government retains for itself, and others acting on its behalf, a paidup nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The U.S. Department of Energy (DOE) will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doepublicaccessplan).
Coupled multimodel simulations often involve high degrees of computationally complex workflows, and achieving consistently accurate solutions is strongly dependent on the choice of spatiotemporal numerical algorithms used to resolve the interacting scales in physical models. Rigorous spatial coupling between components in such systems involves field transformations and communication of data across multiresolution grids while preserving key attributes of interest such as global integrals and local features, which is usually referred to as the process of remapping (also “regridding” or just “interpolation”) (Van Leer, 1979; Dukowicz and Kodis, 1987; Jones, 1999). Such remap procedures are critical in ensuring the stability and accuracy of scientific codes simulating multiphysics problems that typically occur in many different scientific domains. While there have been many highorder, stable interpolators proposed within different disciplines (Zienkiewicz and Zhu, 1992; Grandy, 1999; Jones, 1999; Smith et al., 2000; Garimella et al., 2007; de Boer et al., 2008; Slattery, 2016), to our knowledge, an effort to uniformly compare and measure the properties of these schemes has not yet been attempted. The current work is motivated by a need for an intercomparison of remapping schemes, which led us to standardize several numerical metrics to uniformly compare the properties of these algorithms that are routinely applied to problems in climate and weather system modeling.
Remapping algorithms, in general, compute the spatial interpolation or quasiinterpolation of field data that are defined on a source mesh (Ω_{s}) onto a target (Ω_{t}) component mesh. It is imperative to note that these mesh pairs, Ω_{s} and Ω_{t}, are generally of a different tessellation topology, spatial resolution, and regularity. Over the past several decades, there have been considerable efforts to create new conservative remapping algorithms to tackle various grid combinations and create consistently accurate schemes for solution transfer between component models in multiphysics problems in many fields besides climate. These efforts have resulted in several software libraries and packages developed for this purpose. Examples include schemes developed for fluid–structure interaction (FSI) or heat transfer (such as MpCCI, Joppich and Kürschner, 2006 and preCICE, Bungartz et al., 2016), moving mesh problems with arbitrary Lagrangian–Eulerian (ALE) methods (Dukowicz and Kodis, 1987; Dukowicz and Baumgardner, 2000), and generalpurpose remap software such as MOAB (Tautges and Caceres, 2009; Mahadevan et al., 2020), PANG (Gander and Japhet, 2013), and Data Transfer Kit (DTK) (Slattery et al., 2013).
Remapping packages developed for ESMs, such as the Community Earth System Model (CESM) (Hurrell et al., 2013) and the Energy Exascale Earth System Model (E3SM) (E3SM Project, 2018), include SCRIP (Jones, 1999), the Earth System Modeling Framework (ESMF) Regridder (Hill et al., 2004), TempestRemap (Ullrich and Taylor, 2015; Ullrich et al., 2016), ncremap (Zender, 2008), and the Common Remapping Software (CoR) (Liu et al., 2018) libraries. These remapping packages generate files with mapping weights that are then consumed by couplers such as MCT (Larson et al., 2005), cpl7 (Craig et al., 2012), and OASIS3MCT (Craig et al., 2017), which apply the weights during runtime. Support for generation of mapping weights, and application at runtime is growing (Mahadevan et al., 2020). Note that more recently, a comparative study of various productionready remappers for ESMs using realistic meshes was performed by Valcke et al. (2022) using the specific numerical metrics that have been developed and presented in this paper. The success of this study in better understanding the numerical behavior of production regridders is promising.
Many of the productionready remapping software implementations used in global climate simulations are typically based on first and secondorder conservative meshbased schemes, with additional support for secondorder nonconservative bilinear patch reconstructions (Zienkiewicz and Zhu, 1992; Rider, 2014) or thirdorder bicubic spline interpolations (Hanke et al., 2016) for scalar fields. The projection algorithms implemented in these libraries for the conservative maps are based on the computation of an overlay or a “supermesh” (Jiao and Heath, 2004a; Farrell et al., 2009; Farrell and Maddison, 2011), defined as ${\mathrm{\Omega}}_{S\cap T}={\mathrm{\Omega}}_{\mathrm{s}}\cap {\mathrm{\Omega}}_{\mathrm{t}}$, which is then used to compute consistent and conservative linear operators through L^{2} minimization for transferring field information (Ullrich and Taylor, 2015; Ullrich et al., 2016).
Remapping also occurs in other parts of a climate model, such as within the geophysical fluid dynamics solver of a component. Remapping strategies for tracer transport advection schemes such as the SemiLagrangian Inherently Conserving and Efficient (SLICE) scheme (Zerroukat et al., 2005), Conservative and Monotone Cascade Remap on Sphere for CubedSphere (CS) to regular latitude–longitude (RLL) grids (Lauritzen and Nair, 2008), geometrically exact conservative remapping approaches (Ullrich et al., 2009), and some other highorder massconserving schemes (Norman and Nair, 2008; Erath et al., 2013) have also been devised. While some of the previous work, including conservative semiLagrangian multitracer transport schemes (CSLAM) (Lauritzen et al., 2010) with secondorder variants (Zerroukat et al., 2006), can be easily extended to arbitrary mesh topology combinations, many of the other approaches are tied to specific cases and not used in the intercomponent coupling, which falls out of the scope of this study.
While traditional meshbased schemes requiring computation of overlay–exchange grids have the advantage of being inherently conservative (Balaji et al., 2006), they tend to be computationally expensive (Jansen et al., 1992) and achieving highorder accuracy for scalar field data can sometimes pose difficulties. Numerically suboptimal approximations, such as nodal expansions or nearest neighbor values, can be used to compute interpolations quite efficiently if the field being remapped is smooth and does not mandate conservation prescriptions. Bilinear and bicubic interpolations can also prove to be alternatives here. However, consistently ensuring highorder accuracy under conservation constraints is often necessary for many realworld fields, such as the global heat flux that is exchanged between the atmosphere and ocean components in climate systems.
More recently, new meshbased (Li et al., 2020) and fully meshless (Trask and Kuberry, 2020) remapping schemes that do not require overlap meshes have been developed for other problems in order to mitigate the computational complexity of computing exact intersections between unstructured meshes. These algorithms can provide highorder, 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 realworld simulations. Additionally, while the computational cost is critically important for production runs, our intercomparison study specifically compares only the numerical performance of each algorithm under varying mesh topologies and field regularity, which are closely representative of those from a climate model such as E3SM. The presented protocol hence provides a systematic way to test and compare all existing and new remapping algorithms being developed for Earth system modeling.
Organization
The current study aims to better understand and document the key properties of the remapping schemes through the use of mesh, field, and scale resolutionindependent metrics definitions. This intercomparison paper is organized as follows.

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

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

A sample workflow for computing and comparing various numerical metrics in a unified and unbiased fashion is featured in Sect. 4. Some implementation details on the opensource Pythonbased infrastructure are also provided.

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

Potential future research directions to extend the analysis presented in this work are described in Sect. 6, followed by a summary of the intercomparison study in Sect. 7.
In general, for coupled simulations, we need to transfer a field ${\mathit{\psi}}_{\mathrm{s}}\in {\mathbf{R}}^{{n}_{\mathrm{s}}}$ defined on the source grid, Ω_{s}, to a quantity ${\mathit{\psi}}_{\mathrm{t}}\in {\mathbf{R}}^{{n}_{\mathrm{t}}}$ on the target grid, Ω_{t}. A remap is hence defined as an operator $R:{\mathbf{R}}^{{n}_{\mathrm{s}}}\mapsto {\mathbf{R}}^{{n}_{\mathrm{t}}}$ that transfers ψ_{s} to a trial field ${\stackrel{\mathrm{\u0303}}{\mathit{\psi}}}_{\mathrm{t}}\in {\mathbf{R}}^{{n}_{\mathrm{t}}}$. As noted above, we generally wish to preserve a property defined as an operator $P:{\mathbf{R}}^{{n}_{\mathrm{t}}}\to {\mathbf{R}}^{m}$ that describes discrete physical invariants and/or inequalities relevant to ψ_{t}. However, to keep the current work focused we restrict attention to definitions that are applicable only to scalar fields.
Finally, the accuracy or order of the remap typically depends on how the field ψ is reconstructed from a set of discrete values. On a mesh, this could be a Taylor series expansion around a cell centroid (Jones, 1999) or a general finiteelement (FE) type of expansion using nodal basis functions.
2.1 Related work
One of the simplest data transfer methods is to use piecewise interpolation functionals. This approach is particularly convenient if ψ_{s} defined on Ω_{s} has an associated function space, which can be directly utilized for evaluating the interpolant approximation. Such consistent interpolation techniques utilizing the underlying basis functions for field descriptions give rise to standard secondorder linear (on simplicial elements) or bilinear (with quadrilateral elements) interpolations (Hill et al., 2004). In finitevolume methods or for interpolating from station data (“scattered data” or “location streams”), the nearestpoint interpolation is sometimes used, which is at best firstorder accurate.
More generally, the remap methods designed for scattered data and celltocell 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 highorder interpolators use spline quasiinterpolants (de Boor, 1990), bicubic splines (Hanke et al., 2016; Craig et al., 2017), the standard radialbasis function spaces (Flyer and Wright, 2007; Bungartz et al., 2016), the moving least squares (MLS) method (Lancaster and Salkauskas, 1981), and MLS variants such as the modified MLS (MMLS) (Joldes et al., 2015; Slattery, 2016), which originate from highorder reconstruction methods. Recent extensions to MLS for producing efficient highorder remap involve locally reconstructing the manifold geometry from a point set representation and then generating a compact stencil in the local coordinate chart (Liang and Zhao, 2013; Suchde and Kuhnert, 2019; Trask and Kuberry, 2020; Gross et al., 2020).
More recently, Li et al. (2020) proposed a highorder remap technique for piecewise smooth functions on surfaces, known as WLSENO remap (WLSENOR), which is based on the continuous moving frame (CMF) for smooth functions (Jiao and Wang, 2012) as well as an ENOstyle 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 𝒞^{0} normal field instead of global stencils and 𝒞^{∞} moving frames.
Typically, a remap method can and should be independent of the discretization methods used in physics models. In this context, the method may be nonconservative or conservative with respect to certain properties of the reconstruction. A major disadvantage of the quasiinterpolation techniques is the lack of inherent constraints for conserving energy during the remap process; see, e.g., de Boer et al. (2008) and Jiao and Heath (2004a). Note that the highorder quasiinterpolants (with degree of basis expansion p>1) tend to be significantly less dissipative than loworder methods ($p\in [\mathrm{0},\mathrm{1}]$); see, e.g., Slattery (2016). However, such remap methods are not guaranteed to produce conservative or monotone remapped solution fields.
A common approach to overcome this issue is to enforce conservation so that the integrals over the source and target meshes are equal. There are several examples of first and secondorder conservative remap schemes (Chesshire and Henshaw, 1994; Grandy, 1999; Gander and Japhet, 2013; Jones, 1999; Ullrich and Taylor, 2015) that rely on commonrefinementbased 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 algorithms for efficiency and robustness (Jiao and Heath, 2004b, c). Although these first and secondorder 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 highorder remapping methods, discontinuities in the function defined on Ω_{s} can lead to overshoots or undershoots when evaluated on Ω_{t}. The discontinuities may be in the function values (also called 𝒞^{0} discontinuities), which tend to lead to 𝒪(1) oscillations (or “ringing”) that do not vanish under mesh refinement, analogous to Gibbs phenomena (Gottlieb and Shu, 1997; Fornberg and Flyer, 2007). Additionally, any discontinuities in the derivatives (𝒞^{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 nonphysical variations that can influence the numerical stability of the coupled global nonlinear multiphysics system.
The deficiency in linear mapping approaches arises from the fact that they are only dependent on Ω_{s} and Ω_{t} but completely independent of the solution field that is being projected between the meshes. As a consequence of Godunov's theorem (Godunov, 1959) extended to linear remapping workflows with monotonicity constraints, there may be a restriction on the optimal achievable order of accuracy as shown by Van Leer (1979) while still preserving global solution bounds during the projection step. Hence, the properties of the linear maps can often be enhanced by using a procedure that is nonlinear (depending on the projected field variations) to recover property preservation for high gradient fields.
In this vein, the techniques for resolving Gibbs phenomena can be classified as filtering and mollification; see Jerri (2013). The filtering techniques often rely on postprocessing, such as cropping and property redistribution, to ensure local conservation in a neighborhood. A recent variation of the mass borrowing approach (Royer, 1986) uses a limiter as a filter during the remapping process (Bradley et al., 2019) in order to impose local bounds preservation for linear map applications such that monotonicity can be recovered even when the underlying remapper does not provide it. This “ClipAndAssuredSum” (CAAS) postprocessing filter can be useful to avoid spurious numerical oscillations due to resolution disparity or strong gradients in the underlying solution. Applying CAAS filters to quasiinterpolatory linear maps can hence produce strictly bounded reconstructions on Ω_{t}, while providing property preservation, as a viable option to achieve better remapped solutions in production climate simulations. In contrast, mollification adapts the kernels (i.e., basis functions) near discontinuities. In the past, discontinuity detecting and a posteriori stabilization procedures have been used with specific mesh discretizations (Blanchard and Loubere, 2016) to choose optimal orders of reconstruction adaptively in order to ensure better behavior for polygonal meshes. Other methods, such as WLSENOR, 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 postprocessing 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 optimizationbased 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 “divideandconquer” alternative (Bochev et al., 2014) to direct property preserving methods, which separates accuracy considerations from property preservation. In so doing, OBR helps to avoid interdependencies between mesh quality, accuracy, and property preservation that force many monotone reconstruction methods to make tradeoffs between the latter two (Berger et al., 2005). Extensions of such schemes for remapping climate data are an unexplored topic and may provide avenues for future research.
Additionally, mimetic schemes that use compatible function spaces (Thuburn et al., 2009) depending on the fields being transferred between component models have proven to be extremely accurate for remapping scalar and vector fields on Arakawa C/D grids (Pletzer and Hayek, 2019). Potential extensions for arbitrary mesh topologies to yield compatible, conservative remaps for fluxes in climate components are possible (Ringler et al., 2010). However, to our knowledge, general theory and implementations for remapping on arbitrary meshes are currently unavailable, which may restrict the usage of such schemes for production climate simulations.
2.2 Weighted least squares approximations in remapping
A common theme across all of the remapping methods described in this paper is that they utilize some variants of the least squares approximations (also called quasiinterpolation) in their computational kernels. To illustrate the idea, let us consider a function $f\left(\mathit{u}\right):\mathrm{\Omega}\to \mathbb{R}$ at a given point ${\mathit{u}}_{\mathrm{0}}={\left[\mathrm{0},\mathrm{0}\right]}^{T}$, such as a quadrature point in a cell on the target mesh. Let us first assume a domain Ω⊂ℝ^{2} for simplicity, where $\mathit{u}\equiv [u,v]$. Let f be 𝒞^{p−1} continuous for some degree p≥0, and then f(u) can be approximated to p+1storder accuracy about u_{0} using a degreep twodimensional Taylor polynomial as
where ${c}_{jk}={\displaystyle \frac{\mathrm{1}}{j\mathrm{!}k\mathrm{!}}}{\displaystyle \frac{{\partial}^{j+k}}{\partial {u}^{j}\partial {v}^{k}}}f\left({\mathit{u}}_{\mathrm{0}}\right)$, and h is a measure of the radius of the local neighborhood. In celltocell transfer, some integrals of f are typically known on the source mesh. Given m cells on the source mesh in a neighborhood of a point u_{0} on the target mesh, let ϕ_{i}(u) be the test function (such as the Heaviside functions) associated with the ith cell. Let b_{i} denote these known integral values. We then obtain a system of m equations with $n=(p+\mathrm{1})(p+\mathrm{2})/\mathrm{2}$ unknown coefficients c_{jk},
for $\mathrm{1}\le i\le m$. Note that one can also use a bidegreep Taylor polynomial by letting $\mathrm{0}\le j\le p$ and $\mathrm{0}\le k\le p$ in Eq. (1), which would lead to $n=(p+\mathrm{1})(p+\mathrm{1})$ unknowns. The resulting approximate Taylor polynomial can then be used, for example, to evaluate (or quasiinterpolate) f at u_{0} to p+1st order accuracy or even to higher order due to superconvergence. The same procedure can also be applied to construct a trivariate quasiinterpolation by replacing u and u^{j}v^{k} with x and x^{j}y^{k}z^{ℓ}, respectively. The m equations in Eq. (2) can be written in matrix form as Ax≈b, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ is known as a generalized Vandermonde matrix, x∈ℝ^{n} is composed of c_{jk} in Eq. (1), and b∈ℝ^{m} is composed of the known integrals b_{i} on the source mesh. In general, the generalized Vandermonde system (2) is rectangular. It can be solved by minimizing a weighted norm of the residual vector $\mathit{r}=\mathit{b}\mathbf{A}\mathit{x}$, i.e.,
where W is a diagonal matrix containing the weight for each row in A. This formulation is known as weighted least squares or WLS for short (Golub and Van Loan, 2013). When m=n and A and W are nonsingular, W does not affect the solution. More generally, when m≠n, different W leads to different solutions by changing the relative importance of certain sample points. Different remapping algorithms and reconstruction schemes use specific weighting strategies to achieve optimal solutions to the weighted least squares problem in Eq. (3). Note that the generalized Vandermonde matrix A may be illconditioned as the polynomial degree p used for the reconstruction grows. A preferred approach to address this potential illconditioning is to use a rankrevealing QR factorization (Golub and Van Loan, 2013; Li et al., 2020).
2.3 Algorithms for Earth system models
Among the standard conservative remapping and highorder reconstruction strategies introduced in Sect. 2.1 for climate problems, we selected four specific algorithms to explore in detail. The motivation for choosing these algorithmic implementations was driven by their high usability, including the use in current ESMs, and the rigorous underlying numerics that can be verified and validated consistently for a large suite of test problems. We categorize these algorithms below into three broad groups based on whether the algorithms require overlay meshes and whether mesh data structures are utilized to compute the solution reconstruction.
 I.
Overlaymeshbased 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 higherorder L^{2} projection support
Both ESMF Regrid and TempestRemap provide the remapping capability for scalar fields defined on Ω_{s} based on a weighted residual formulation given in Eq. (3). With a source field ψ_{s} defined on Ω_{s}, the formulation computes ψ_{t} on Ω_{t} by solving the following problem:
$$\begin{array}{}\text{(4)}& \underset{\mathrm{\Omega}}{\int}{\mathit{\psi}}_{\mathrm{t}}{\mathit{\varphi}}_{i}\phantom{\rule{0.33em}{0ex}}\mathrm{d}V=\underset{\mathrm{\Omega}}{\int}{\mathit{\psi}}_{\mathrm{s}}{\mathit{\varphi}}_{i}\phantom{\rule{0.33em}{0ex}}\mathrm{d}V,\end{array}$$where the ϕ_{i} represents suitable weight functions. In a commonrefinementbased transfer (Jiao and Heath, 2004a), the ϕ_{i} represents the basis functions ψ_{t,i} on Ω_{t}, which leads to a Galerkin projection method. Such a Galerkin method can achieve conservation globally by solving the weighted residual minimization in Eq. (4). More details on the specific methodologies used for reconstruction in these implementations are provided in later sections.
Note that the SCRIP library (Jones, 1999) also provides both first and secondorder conservative remapping schemes, but it was not included in the current study since its capabilities are similar to that of ESMF and TempestRemap implementations, and those two are used in current production versions of CESM and E3SM.
 a.
 II.
Meshless remappers
Reconstruction methods that do not directly utilize the topological information about the underlying mesh layout are meshless methods. In our study, we will consider the generalized moving least squares (GMLS) method, with global conservation, monotonicity constraints, and local bounds preservation provided by CAAS as a postprocessing 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.
Nonoverlay meshbased hybrid remappers
The final category includes the meshbased remappers that do not require computation of an intersection mesh between Ω_{s} and Ω_{t}. The weighted least squares essentially nonoscillatory remap method (WLSENOR) utilizes a discontinuitycapturing, highorder technique for piecewise smooth functions to produce optimal field transformations between meshes by minimizing the residual in Eq. (3).
Other opportunities for comparison in this category include using reconstructed climate data with the conservative bilinear algorithm and patch reconstructions in ESMF or bicubic interpolations available in Yet Another Coupler (YAC) (Hanke et al., 2016).
These algorithms and their implementations span a range of remapping techniques, including those currently used in production runs and those that can potentially be used in the future given the availability of opensource software. Further details regarding the numerics and the software tools for each of the schemes are provided in the following subsections.
2.3.1 ESMF
The Earth System Modeling Framework (ESMF, https://earthsystemmodeling.org/, last access: 7 August 2022) (Hill et al., 2004) consists of a suite of software modeling tools that support building Earth system models. Among other features, ESMF exposes several key functionalities to transfer data between component models in weather and climate applications. It offers a variety of data structures for transferring field data between components and libraries for regridding, time advancement, and other common modeling functions. The advanced regridding algorithms provided by ESMF implement standard bilinear interpolation, higherorder patch recovery, first and secondorder conservative projections, and several variations of nonconservative nearest neighbor interpolants. ESMF can produce maps that are “offline” in the sense that they can be precomputed, stored, and then applied as a linear operators to arbitrary fields defined on Ω_{s} to compute the projection on Ω_{t}. Fully online remapping is also possible with ESMF.
In the current intercomparison study, we utilized the commandline applications installed with ESMF (version 8.1) to generate and apply interpolation weights from the command line using NetCDF files. Using this ESMF_RegridWeightGen application, interpolation weights can be generated for various grid combinations and applied to source field data to compute projections between component models. More specifically, we select only the first and secondorder conservative projection methods invoked with “conserve” and “conserve2nd” commandline options, respectively, which are implemented in ESMF for scalar fields and fluxes. These options are routinely used in production climate models and hence provide valuable baselines on the current state of remapping algorithms in our comparison study. For more details regarding the numerics and implementation of the algorithms in ESMF, we refer readers to Collins et al. (2005), Jones (1999), and Kritsikis et al. (2017).
2.3.2 TempestRemap
TempestRemap (https://github.com/ClimateGlobalChange/tempestremap, last access: 7 August 2022) is a software library to generate conservative, consistent, and monotone remapping weights of arbitrarily highorder accuracy (with degree p>1) between meshes on the sphere. Similar to the ESMF tools, TempestRemap generates offline maps that are consumed by ESMs. TempestRemap supports the computation of remapping weights between any combination of finitevolume (FV) and/or spectral element (SE) (both continuous and discontinuous spaces) descriptions with conservative prescriptions. For purposes of the current paper, we describe only the highorder FV–FV algorithms below.
In TempestRemap, the procedure used to generate remapping weights for FV discretizations consists of two primary operations (Ullrich and Taylor, 2015; Ullrich et al., 2016). First, local polynomial reconstructions are defined over the source mesh with some adjustments made for spherical geometry (Jalali and Gooch, 2013). The coefficients of these local reconstructions are computed according to a weighted 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 ${\mathrm{\Omega}}_{\mathrm{s}\cap \mathrm{t}}={\mathrm{\Omega}}_{\mathrm{s}}\cap {\mathrm{\Omega}}_{\mathrm{t}}$. Note that when computing highorder linear maps, if a sufficiently large patch on the source map is not used during the reconstruction process, the condition number of the generalized Vandermonde matrix A in Eq. (2) can be very high, leading to numerical roundoff errors and degradation in the overall accuracy of the remap operator to firstorder accuracy in the neighborhood. Further details on the analysis are provided in Ullrich and Taylor (2015).
A common approach for the integration operator has been to construct a potential function with divergence equal to the scalar field being integrated and then use the divergence theorem to transform integrals over mesh faces into line integrals around the face (Dukowicz and Kodis, 1987; Jones, 1999; Ullrich et al., 2009). While this technique has been used to define conservative remapping operators, a drawback is that it requires an analytical expression for the potential function, which can be difficult in general. An alternative method was proposed in Erath et al. (2013), which, although avoiding some of the difficulties of the line integral approach, lacks consistency. In TempestRemap, a quadraturebased integration operator is used to generate an initial set of remapping coefficients, which is guaranteed to be consistent (exactly remaps the unit field) but may not be conservative. To obtain conservation, the coefficients of the matrix are projected linearly into the space of maps that are both conservative and consistent.
2.3.3 GMLS
The generalized moving least squares (GMLS) method extends the moving least squares (MLS) technique from approximation of point values to approximation of arbitrary linear functionals (Nayroles et al., 1992; Breitkopf et al., 2000; Wendland, 2004; Mirzaei et al., 2012). Similar to MLS, a distancebased weighting kernel is used to favor the source data sites closer to the point of query or reconstruction. GMLS features a broad choice of sampling functionals and target operators, and hence the term “generalized” is in its name. The application of nonlinear sampling functionals and target operators for GMLS is possible (Gross et al., 2020) but is certainly not common. In fact, most applications of GMLS use linear sampling functionals and target operators (gradient, curl, divergence, integral average along an edge or over a cell, etc.), for which there is theory on approximation and wellposedness. When the sampling functionals and the target operator are simply pointwise evaluations, GMLS reduces to the traditional MLS method.
Highorder accurate approximations cannot be achieved with traditional MLS schemes using function spaces described by Raviart–Thomas (ℋ(div)) and Nedelec basis (ℋ(curl)) or even cellaveraged quantities that are common in FV codes. However, through careful selection of the sampling functionals and target operators, consistent approximations of fields embedded in these various nonstandard spaces are possible using the GMLS approach. This flexibility of the GMLS method greatly increases the available types of input data that can be handled to produce highorder reconstructions of fields between Ω_{s} and Ω_{t}. When a GMLS problem is solved, a set of coefficients corresponding to elements of a basis (traditionally polynomial) is computed. The target operation is then applied to each member of the basis used for approximation, after which a dot product is made between the originally computed coefficients and the evaluation of the target operation acting on the basis. While GMLS permits a wide choice of target operators, we choose a target operation which is the cell average of each member of the approximating basis on the target grid.
Traditionally, GMLS uses a basis that is defined as a function of the spatial dimension from which a point cloud is sampled. However, in this work, reconstruction of functions sampled on a manifold permits generating a compact stencil in a local coordinate chart, which is one dimension smaller (Liang and Zhao, 2013; Suchde and Kuhnert, 2019; Trask and Kuberry, 2020; Gross et al., 2020). The savings in net computational floatingpoint operations (flops) is a factor of p^{3} in R^{2} 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 discretizationdependent and generally vanishes with refinement. However, such weak conservation notions may not be deemed satisfactory for climate modeling applications for which exact global conservation has a history of being demanded and valued. In such cases, GMLS remap should be followed by a postprocess filter to restore global conservation to machine precision. In this paper, we use either the GMLS or GMLSCAAS notations to indicate whether the CAAS routine has been used as a postprocessing filter after each remap step in order to restore global conservation or global bounds, along with an attempt to improve local property preservation.
Similar to the overlaymeshbased, ESMF, and TempestRemap offline remappers described in Sect. 2.3.1 and 2.3.2, respectively, the solution of the GMLS problem produces a stencil or a linear mapping that can be computed a priori to a simulation. This enables storage of the relevant parts of the remap operator as a sparse matrix. Therefore, the use of GMLS without postprocess filtering can be thought of as an offline remapper, and GMLSCAAS can combine the offline and online processes to achieve more favorable properties. Note that the CAAS filter is only one of several choices available for the online postprocessing filter. Alternative strategies such as the nonlinear OBR (Bochev et al., 2014) for feature detection and minimization of local dissipation may be possible in conjunction with the GMLS workflow. These enhanced GMLS remapping strategies may be considered in future studies.
There are several key motivations for using GMLS to perform field remapping for ESMs. These include mesh topology independence, as well as flexibility in the choice of the sampling functional and the target operator, thereby enabling remap of nontraditional degrees of freedom that may be defined on the vertices, edges, or faces of Ω_{s}. Additionally, the embarrassingly parallel nature of the dense linear systems that require inversions during reconstruction yields a high flopstocommunication ratio. This feature makes GMLS more suited for nextgeneration platforms. Given the computationally intensive nature of the GMLS algorithm, it is best implemented in libraries focused on parallel performance portability, computational efficiency, and significant compiler optimizations. The Compadre Toolkit (https://github.com/SNLComputation/compadre, last access: 7 August 2022) (Kuberry et al., 2019) is built on Kokkos (Edwards et al., 2014) and Kokkos kernels such that a single version of the code is written to be performant on both multicore CPU and hybrid GPU architectures. Compadre is built to assemble and solve large batches of GMLS problems in parallel, thereby leveraging batched QR solvers with pivoting for the parallel solution of many small dense linear systems.
2.3.4 WLSENOR
WLSENOR, or weightedleastsquaresbased essentially nonoscillatory remap (Li et al., 2020), is based on the WLSENO reconstruction technique proposed in Liu and Jiao (2016). Originally developed for solving hyperbolic partial differential equations (PDEs), WLSENO detects discontinuities and then reduces or eliminates the potential Gibbs phenomena in the solutions of the PDEs by adapting the weighting schemes in WLS based on the local features in the solution. WLSENOR adapted WLSENO to remap data between meshes, and in the process, it enhanced the treatment of discontinuities to resolve not only the severe oscillations at 𝒞^{0} discontinuities (i.e., jumps in function values), but also the accumulated effect of mild oscillations due to 𝒞^{1} discontinuities (i.e., jumps in the derivatives).
Unlike the preceding remapping techniques, WLSENOR is a nonoverlay meshbased technique in that it uses the mesh for computing numerical integration, but it does not require an overlay mesh (although it has the option to use the overlay mesh if available). More specifically, WLSENOR utilizes adaptive quadrature rules with p refinement in smooth regions and h refinement near discontinuities in its numerical integration to achieve highorder accuracy (with degree p>1) and (nearly local) conservation. More specifically, WLSENOR is composed of three major components. The first component is a WLSbased algorithm for smooth functions, as we described in Sect. 2.2. In this context, the weighting scheme in WLSENOR is based on a positivedefinite radial function due to Buhmann (Buhmann, 2001). As shown in Li et al. (2020), this weighting scheme encourages statistical error cancelations and enables better superconvergence (i.e., higher than (p+1)storder convergence) with evendegree p for nodetonode transfer of smooth functions. In this intercomparison work, we use an extension of the work in Li et al. (2020) to cellcentered data. Mathematically, this extension essentially uses the step functions (also called Heaviside functions) over the cells on the source mesh as the test functions in Eq. (2) compared to the Dirac delta function at the nodes on the source mesh as test functions in Li et al. (2020). Note that in this work, we apply WLS on a sphere Ω=S, which is topologically a 2D object embedded in ℝ^{3}. One could construct a WLS in ℝ^{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 firstorder estimation. 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 quasiinterpolation. In terms of implementation, WLSENOR constructs a matrixbased transfer operator for smooth functions, which maps the cellaveraged values from the source mesh to the target mesh.
The second component in WLSENOR is the detection and resolution of discontinuities. In particular, WLSENOR can detect 𝒞^{0} and 𝒞^{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 WLSENOR is based on an asymptotic analysis of the Gibbs phenomena near 𝒞^{0} and 𝒞^{1} discontinuities as described in Li et al. (2020). The original detector in Li et al. (2020), however, was for nodetonode transfer. For celltocell transfer in this work, we simply apply the detector to the dual mesh and treat cellaveraged values as approximations of cellcenter values. This approximation is secondorder accurate, which is sufficient in detecting discontinuities. After identifying the discontinuities, WLSENOR uses a solutionbased weighting scheme that effectively leads to (nearly) onesided quasiinterpolation in discontinuous regions. This solutionbased weighting scheme causes WLSENOR to use a different set of basis functions to overcome the Gibbs phenomena, and hence it can be classified as a mollification technique (see Sect. 2.1). In terms of implementation, the discontinuity detector on the source mesh is implemented as a matrixbased operator; after the discontinuities are identified, a new local solutionbased 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 WLSENOR is an adaptive quadrature technique, which is enabled when the target mesh is significantly coarser than the source mesh. In this case, simply sampling the function at the quadrature points of a target cell may miss some important local features on the source mesh, especially near discontinuities. To overcome this loss of information, one could use the overlay mesh as in TempestRemap. This approach may be ideal in terms of accuracy and conservation, but it introduces complications when the elements have curved edges. Although WLSENOR implementation supports this option, in this comparative work, we use a nonoverlaybased version of WLSENOR that utilizes h and p refinement of the quadrature rules in discontinuous and smooth regions, respectively. In particular, near the detected discontinuities, WLSENOR subdivides the cells (i.e., hrefinement) on the target mesh to match the local resolution on the source mesh and then utilizes the quadrature points of the subdivided cells. For smooth regions, we found it more efficient to use the quadrature points of higherdegree quadrature rules (i.e., p refinement) than subdividing the cells. This adaptive quadrature technique not only overcomes the loss of accuracy but also enables WLSENOR to recover global conservation in a fashion that is nearly local to discontinuous regions. To this end, if there is a gain or loss in terms of global conservation, we distribute this global error proportionally to cells that have gained or lost mass locally, correspondingly. This conservation–recovery procedure reduces the pollution of the conservation errors from discontinuous regions to smooth regions. We estimate the local conservation errors using the subdivided cells near discontinuities for accuracy and robustness; for smooth regions, we use a simple comparison of the local extreme values in the local neighborhood for smooth regions for efficiency.
The current implementation of the WLSENOR algorithm uses MATLAB, with which the core components were converted into C using “MATLAB Coder” (version 4.2). An opensource C++ implementation is currently underway and will be released in the future for both nodetonode and celltocell field transfers.
Solution remapping on unstructured meshes is a complex process, and it is critical to satisfy several key properties to ensure that the transferred field data between components do not introduce unbounded and nonphysical error modes. In order to compare different remapping schemes in an unbiased framework, we introduce five primary categories under which the comparison metrics can be grouped.

Sensitivity is algorithmic invariance to underlying component mesh topology.

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

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

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

Dissipation is the minimization of local solution dispersion on repeated backandforth remap transfers between Ω_{s} and Ω_{t}.
Given a continuous field ψ, we use D^{s} and D^{t} to denote the reference discretization of ψ on the source and target grid, respectively. These are generated by directly sampling the analytical fields and by spherical harmonic expansion (see Sect. 4) for the realistic fields. The regridding operator from the source to target mesh is denoted by R; i.e., the regridded field is denoted by RD^{s}[ψ]. The global integral operator is denoted by I^{s} and I^{t} on the source and target grid, respectively. Typically these operators take the form
where J_{j} denotes the area or weight of the jth degree of freedom (i.e., the volume of a finite volume).
3.1 Grid sensitivity
A crucial factor for the success and broader usability of a general remapping algorithm in ESMs is the ability to produce 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 to compare remapping schemes using the cubedsphere (CS), quasiuniform Voronoi (MPAS), regular latitude–longitude (RLL) meshes on both quasiuniform and regionally refined meshes (RRMs). Some sample meshes used in the study are shown in Fig. 1.
3.2 Standard accuracy measures
Accuracy in the remapped solution will be assessed with standard error metrics defined as follows.
where ${\mathit{e}}_{K}=\mathrm{}{\mathbf{RD}}^{\mathrm{s}}\left[{\mathit{\psi}}_{K}\right]{\mathbf{D}}^{\mathrm{t}}\left[{\mathit{\psi}}_{K}\right]\mathrm{}$ is the local remapping error in element K∈Ω_{t}, E represents the discrete error vector in the solution field relative to the reference field sampled on Ω_{t}, and I^{t} is the weighted integral using Eq. (5) on Ω_{t}. In some general sense, the error measure ${\u2225E\u2225}_{{L}_{\mathrm{1}}}$ identifies errors in largescale features, ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ identifies errors in smallscale features, and ${\u2225E\u2225}_{{L}_{\mathrm{\infty}}}$ identifies the largest pointwise error. These are special cases of error measures ${\u2225E\u2225}_{{L}_{p}}$ in L_{p} norm, in which as p tends from 1 to ∞, the norms capture the largescale mean absolute errors and the largest pointwise errors at the two extremes, respectively, for p=1 and p=∞. Related to accuracy, consistency is assessed by applying these metrics to uniformly smooth fields with no C^{0} or C^{1} discontinuities and verifying the asymptotic theoretical rate of convergence under uniform mesh refinement conditions. The asymptotic convergence order of a given remapping method with a degreep reconstruction is in general expected to be 𝒪(h^{p+1}) in ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ and 𝒪(h^{p}) in ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$, where h is the characteristic spatial length of the mesh that can be computed using a simple definition such as $h=\frac{\mathrm{1}}{{N}_{{\mathrm{\Omega}}_{\mathrm{t}}}\left(\mathrm{el}\right)}$, with ${N}_{{\mathrm{\Omega}}_{\mathrm{t}}}\left(\mathrm{el}\right)$ representing the number of elements in the unstructured mesh Ω_{t}. Note that for nonuniform meshes such as RRMs, it is important to take the root mean square (rms) element Jacobian value for the parameter h. Additionally, while the accuracy of the remap projection is strongly dependent on both Ω_{s} and Ω_{t}, in studies presented here, we only utilize uniform refinements in both meshes (without including crossresolution calculations resulting from a fine Ω_{s} and coarse Ω_{t}) to compute the overall convergence order for projecting smooth fields.
Note that in order to eliminate potential aliasing errors, the normalization factors D^{t}[ψ] used in the denominator for definitions of ${\u2225E\u2225}_{{L}_{\mathrm{1}}}$, ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$, and ${\u2225E\u2225}_{{L}_{\mathrm{\infty}}}$ are computed based on the exact sampling (elementaveraged for FV discretization) of the data (reference solution) on Ω_{t}, and not using the projection of the field RD^{s}[ψ]∈Ω_{t}.
3.3 Gradient preservation measures
Preservation of the solution gradients in addition to other critical properties, such as local conservation in the remapping procedure, requires C^{1} continuity in the D^{s}[ψ]. Let ∇D^{s} and ∇RD^{t} be the gradients of the scalar fields on Ω_{s} and Ω_{t}, respectively.
Then, in order to measure accuracy of the solution and its gradient, we introduce two specific global metrics: ${\leftE\right}_{{H}_{\mathrm{1}}}$ seminorm and the ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ norm.
where ${\mathit{e}}_{K}=\mathrm{}{\mathbf{RD}}^{\mathrm{s}}\left[{\mathit{\psi}}_{K}\right]{\mathbf{D}}^{\mathrm{t}}\left[{\mathit{\psi}}_{K}\right]\mathrm{}$ is the local remapping error in element K, $\mathrm{\nabla}{\mathit{e}}_{K}=\mathrm{}\mathrm{\nabla}{\mathbf{RD}}^{\mathrm{s}}\left[{\mathit{\psi}}_{K}\right]\mathrm{\nabla}{\mathbf{D}}^{\mathrm{t}}\left[{\mathit{\psi}}_{K}\right]\mathrm{}$ is the corresponding gradient of the error, and I^{t} is the weighted integral using Eq. (5) on Ω_{t}.
In this study, ∇q for some scalar q, defined as a cell mean value, is generated by finding the convex hull of cells surrounding the current cell and computing the gradient per Barth and Jespersen (1989). We compute these gradients as part of the metrics evaluation after each remapping sequence.
3.4 Global conservation
Global conservation is trivially assessed by evaluating the change in the global integral of the scalar field value on the source mesh and the projected field on the target mesh. We use the following metric to quantify global conservation.
However, we note that this definition for L_{g} is only meaningful when the target domain fully envelops the source domain (which may have gaps or holes in more general cases). In the case of climate modeling, an admissible example for using Eq. (11) would be for remapping heat and moisture fluxes from the land surface to the overlying atmosphere.
3.5 Global extrema preservation
Global extrema preservation can be assessed via the standard G_{min} and G_{max} error metrics (Ullrich and Taylor, 2015):
The error measures G_{min} and G_{max} identify undershoots and overshoots, respectively, by taking on nonzero values ($\left{G}_{\mathrm{min}}\right\le \mathrm{0}$ and $\left{G}_{\mathrm{max}}\right\ge \mathrm{0}$) when there is a departure away from the reference global extreme values. In other words, a nonzero value of the metric indicates changes in global extrema, indicating the presence of a Gibbs phenomenon or instabilities introduced due to the interpolation. Hence, the global extrema metric is particularly useful as it provides indications about the monotonicitypreserving properties of the remapping schemes.
3.6 Local extrema preservation
Local extrema preservation can be assessed using a localized difference; i.e., to what degree does the remapped grid cell value fall within the range of surrounding grid cells sampled on the target grid? This consideration motivates us to define a localized difference in extrema:
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 Eqs. (14) and (15) utilizes a local neighborhood to determine the deviation from reference extrema values. This is sufficient to capture resolution of sharp gradients in the remapped fields under mesh refinement for elementaveraged data. However, the metric contains 𝒪(h) dependence on the mesh resolution and can be applied to 𝒞^{0} or smoother fields, but not when 𝒞^{0} discontinuities are present.
For all remapping algorithms evaluated in this comparison study, we conduct iterative twoway (Ω_{s}→Ω_{t} and Ω_{t}→Ω_{s}) remapping of an initial source field with FV discretization on Ω_{s}. We do so in order to characterize the stability of a scheme and expose any dissipation effects, which would not be possible to ascertain when comparing single, unidirectional field transfers. With this workflow, we seek to quantify the consistency, stability, and convergence of each participating algorithm as measured with the metrics defined in Sect. 3. Here, N_{r} indicates the number of iterative applications of the linear map to compute field transformation on ${\mathrm{\Omega}}_{\mathrm{s}}\to {\mathrm{\Omega}}_{\mathrm{t}}\to {\mathrm{\Omega}}_{\mathrm{s}}$, which will be referred to as remap iterations. While production ESM solvers do not utilize repeated remap transfers at every time step, our approach to use iterative remaps can provide valuable insight into the dissipative effects and longterm temporal behavior of fully coupled climate simulations, in addition to determining the stability of the remap operator without a need for explicit spectral analysis.
4.1 Opensource MIRA implementation
The workflow necessary to evaluate a given remapping method comprises five consecutive steps described below.

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

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 highorder 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.

A second Python module called FieldGenDriver then takes each of the preprocessed mesh files and evaluates scalar fields by sampling from either an analytical function on the sphere or a set of prescribed spherical harmonic (SPH) coefficients, which is described in Sect. 4.2. In this step, a cell average is computed by local quadrature within each mesh element to a given order of accuracy by appropriately choosing the order of the quadrature to resolve the SPH expansion order. This operation is performed on all the input meshes to generate the reference “ground truth” realization of a given field, which is used to accurately compute the metrics defined in Sect. 3.
We emphasize that any existing mesh (such as Yin–Yang Kageyama and Sato, 2004, or cubic–octahedral – Gaussian – reduced grids) can be used in this workflow, instead of the ones generated in step 1, as FieldGenDriver only relies on the existence of element connectivities and adjacency maps to be available for computing cell integrals.

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

The final Python module in the metrics suite, MetricsDriver, can then be invoked on each of the remapped output data (typically stored in NetCDF files) to consistently compute all the remapping metrics defined in Sect. 3. The computed remapping metrics are then stored as commaseparated values (CSVs) for further analysis and intercomparison studies.
The schematic shown in Fig. 2 provides further details on this workflow. Note that in order to evaluate and compare a new remapping implementation such as SCRIP (Jones, 1999), CoR (Liu et al., 2018), or YAC (Hanke et al., 2016), only the fourth and fifth steps in the workflow have to be executed, since the preprocessed input meshes and the sampled reference data have been made publicly available (Mahadevan et al., 2021).
4.2 Scalar test variables on the sphere
In the remapping intercomparison study, we consider five scalar test variables defined on the sphere as reference solutions fields. These fields are chosen such that different aspects of the remapper can be evaluated uniformly. Details about the analytical and realworld fields as well as the sampling methodology used in the Python implementation are provided below.
4.2.1 Idealized analytical fields
Idealized fields used in this study mirror the approach of (Lauritzen and Nair, 2008) and (Ullrich and Taylor, 2015). Namely, we employ three idealized test cases of varying complexity to understand the error measures produced by remapping. The two analytical fields studied are depicted in Fig. 3.
The first analytical field (AnalyticalFun1) is a combination of spherical harmonics functions with frequency wave similar to order 3, given by
where ${Y}_{m}^{l}$ represents the real spherical harmonic functions evaluated through the SHTools package for degree m and polynomial order l.
Following (Jones, 1999), and (Lauritzen and Nair, 2008), the second field (AnalyticalFun2) is a relatively smooth function resembling spherical harmonics of order 2 and azimuthal wavenumber 2, given by
These fields are used to test performance for a smooth, wellresolved field and a slightly highfrequency, 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 closedform functions and evaluate them on the sphere by using highorder quadrature order rules to sample and compute elementaveraged data. This design allows the flexibility to test slightly more complex analytical vortex fields (Ullrich et al., 2009) or any threedimensional realvalued function projected on the sphere with coordinate transformations (Townsend et al., 2016).
4.2.2 Real data fields
We also test the performance of each remap technique by regridding real data fields obtained from freely available composite satellite observations. The fields chosen are total precipitable water (TPW), cloud fraction (CFR), and global topography (TOPO). From a representative dataset, we compute a spherical harmonic (SPH) decomposition in order to determine an analytical approximation of spectral content. We employ these particular fields because we can control their characteristics as functions, i.e., global bounds for topography, positive definiteness for precipitable water, and continuity for cloud fraction. As such, these fields present distinct challenges to the remapping methods. For convenience, we use the SHTools (Wieczorech and Meschede, 2018) package through its Python interface (PySHTools), which facilitates the computation of spectra based on spherical harmonic bases and reconstruction of fields thereof.
Given the 1D averaged spatial amplitude spectrum for each set of composite satellite data as shown in Fig. 4 with a corresponding linear fit, we then produce controlled randomized realizations for each field on any unstructured mesh and resolution, including regionally refined grids. The randomization is applied to the coefficients of the expansion at each degree, extracted from the linear fit functions in Fig. 4, and is entirely reproducible given an integer seed to a pseudorandom number generator (set to 384 in our study). The code for this is available in the FieldGenDriver (Guerra et al., 2021). The original data range from 0.25 to 0.1^{∘} of resolution and SPH reconstructions are smoothly varying up to the number of coefficients employed. Note that while the order of SPH reconstruction can be specified arbitrarily high (between 1 and 512 modes) to get a better resolution of the fields, the computation of elementaveraged sampling representative of FV discretization needs to utilize a sufficiently highdegree quadrature rule such that the SPH expansions are exactly integrated on the sphere.
 Total precipitable water (TPW).

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

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

Global topography data are taken from the ETOPO1 Global Relief Model (https://www.ngdc.noaa.gov/mgg/global/global.html, last access: 7 August 2022) (NOAA National Geophysical Data Center, 2009; Amante and Eakins, 2009). This field has a global minimum and maximum for all realizations of −10994.0 and 8848.0 m but is otherwise smoothly varying. A randomized reconstruction of TOPO is shown in Fig. 5 (bottom).
The raw satellite data snapshot of TPW, CFR, and TOPO fields that are used in the current study have been made available separately (Guerra and Mahadevan, 2021) in order to make the workflow reproducible.
Comparing different remapping algorithms under a unified infrastructure for test problems and metrics collection is a nontrivial task. The metrics defined in this study and the implementation of the various field samplings on arbitrary unstructured meshes have provided large output datasets to analyze the key properties of the remappers under consideration.
Specifically, for uniformly refined experiments, a series of different mesh types (CS, MPAS, RLL) ${N}_{\mathrm{type}}^{\mathrm{uni}}=\mathrm{3}$ with five different refinement levels (${N}_{\mathrm{ref}}^{\mathrm{uni}}=\mathrm{5}$) was generated. Similarly, for the regionally refined experiments, two different mesh types (${N}_{\mathrm{type}}^{\mathrm{rrm}}=\mathrm{2}$) using CS and MPAS with three refinement levels (${N}_{\mathrm{ref}}^{\mathrm{rrm}}=\mathrm{3}$) around the continental US were used. Table 1 provides the details on the number of elements and nodes in the various meshes utilized in the current study.
Using the field definitions (N_{fields}=5) introduced in Sect. 4.2, consisting of two smooth analytical fields and three real fields expanded with spherical harmonics, sampling was then evaluated on all input meshes (${N}_{\mathrm{type}}^{\mathrm{uni}}{N}_{\mathrm{ref}}^{\mathrm{uni}}+{N}_{\mathrm{type}}^{\mathrm{rrm}}{N}_{\mathrm{ref}}^{\mathrm{rrm}}=\mathrm{15}+\mathrm{6}=\mathrm{21}$) to serve as reference solutions. With these input meshes (Mahadevan et al., 2021), the remapping algorithms were applied following the workflow in Fig. 2 for combinations of CSMPAS (both uniform and RRM), MPASRLL, and RLLCS meshes of varying resolutions to evaluate the metrics data.
The volume of consolidated output metrics data is enormous from this experiment, since 1000 remapping iterations were performed on $N={N}_{\mathrm{type}}^{\mathrm{uni}}{C}_{\mathrm{2}}[{N}_{\mathrm{ref}}^{\mathrm{uni}}{]}^{\mathrm{2}}{N}_{\mathrm{fields}}=\mathrm{375}$ global, uniformly refined mesh cases and $N={N}_{\mathrm{type}}^{\mathrm{rrm}}{C}_{\mathrm{2}}[{N}_{\mathrm{ref}}^{\mathrm{rrm}}{]}^{\mathrm{2}}{N}_{\mathrm{fields}}=\mathrm{45}$ RRM cases, with each of the four remapping algorithms using various degrees of reconstructions p for FV–FV field transfers. Measuring more than 15 different remapping metrics for each of these cases has provided extensively detailed results to compare the algorithmic implementations in an unbiased fashion. Hence, unless explicitly noted, only significantly unique results are presented and discussed below in the following subsections. We direct readers to the IPython notebooks available in Mahadevan et al. (2021), which can be used to generate the comparison for any combination of mesh types, source resolution, target resolution, field variable, and remapping schemes. Note that all the metrics data collected during the analysis are also stored in the same repository for reproducibility.
Detailed results from the intercomparison study and discussion on the implication of each metric to the remapping scheme are presented next.
5.1 Consistency
The consistency of the highorder remap algorithm implementations can be verified by remapping smooth functions and calculating the spatial convergence order of the resultant approximations on the target mesh after repeated remaps. Using the sampled analytical functions described by AnalyticalFun1 and AnalyticalFun2, a verification study was conducted using the standard error norms and gradient error metrics data for the various schemes. For smooth solution profiles, the theoretically expected convergence rates for a consistent remapping method are 𝒪(h^{p+1}) in ${\u2225E\u2225}_{{L}_{\mathrm{1}}}$, ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$, and ${\u2225E\u2225}_{{L}_{\mathrm{\infty}}}$ global error norms and 𝒪(h^{p}) in ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ and ${\leftE\right}_{{H}_{\mathrm{1}}}$ global gradient error norms. The convergence rates measured in the current studies use the definition of h provided in Sect. 3.2 by using the N_{Ω}(el) data described in Table 1 for different refined meshes.
In general, from these studies we observed that convergence rates for both low and highdegree approximations up to p=2 of the analytically smooth fields show good agreement with the theoretically expected accuracy convergence rates. However, for highorder meshbased remaps, achieving a convergence rate higher than 𝒪(h^{3}) appears to be a limitation with the methods involved in this study, while the hybrid and meshless schemes show consistent recovery of the approximation orders for all degrees tested. The following subsections provide detailed results and discussions for each of the remapping algorithms.
5.1.1 ESMF
The conservative schemes implemented in the ESMF package (Hill et al., 2004) have been thoroughly verified and are routinely used to generate the linear maps for solution transfer between components in E3SM. The secondorder conservative projection algorithm that was originally introduced for ALE computations (Dukowicz and Kodis, 1987), and later applied to spherical meshes (Jones, 1999), has been implemented in ESMF with an appropriate linear gradient reconstruction (Kritsikis et al., 2017). We measure the convergence rates for both the firstorder (conserve) and secondorder (conserve2nd) conservative schemes and present the results observed in Table 2. The ESMF firstorder conservative scheme yields expected rates of 𝒪(h) asymptotically. However, the secondorder scheme (conserve2nd) shows degraded convergence rates as confirmed by the global error and gradient norms. This convergence result is unexpected and contrary to the analysis of the secondorder, piecewise linear finitevolume reconstruction procedure presented by Kritsikis et al. (2017), which is implemented in ESMF.
Note that we have presented the computed convergence rates from the analysis as is, and qualitatively speaking, the values in Table 2 that are approximately equal to 1.0 indicate 𝒪(h) and approximately equal to 0.0 indicate 𝒪(1) behavior.
In order to better understand the relative accuracy of the first and secondorder conservative remapping implementations in ESMF, we performed a comparative analysis for all grid combinations using standard global error norms. These results are shown in Fig. 6, where the ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ and ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ error profiles for both the conservative schemes are compared as a function of remap iterations (N_{r}) for the r4 resolution meshes. The computed error metrics from both the schemes indicate that, even though the convergence rates are similar, the conserve2nd option in ESMF produces a significantly better approximation for the remapped field, irrespective of the mesh resolution or field characteristics in the tested samples. It is also evident that the ESMF algorithms are more accurate for the RLLCS mesh types in all global error norms, indicating sensitivity to the underlying element types in the mesh. The significantly better error approximations of the conserve2nd conservative scheme emphasize that it should be preferred over the conserve scheme when possible, irrespective of the underlying mesh types involved in remap.
5.1.2 TempestRemap
The conservative highorder linear maps computed by TempestRemap, as shown in Table 3, produce higherorder expected theoretical convergence rates in comparison to ESMF for the smooth analytical solution fields. The convergence results presented here have been generated with a bidegreep basis reconstruction using a rectangular truncation strategy that computes source patches based on edgeadjacency graphs.
Note that even for smooth solutions, the convergence rates observed for p>2 degrade during the conservative reconstruction for FV–FV maps. This degradation is evident in the rate reduction in all the global error norms for p>2, which indicates that a larger patch size may be necessary to accurately and consistently recover the smooth field after remap from Ω_{s} to Ω_{t}. The rate reduction is especially noticeable in the presence of pole singularities like that of RLL meshes, where the ${\u2225E\u2225}_{{L}_{\mathrm{\infty}}}$ shows further reduction in convergence rates. The failure to achieve 𝒪(h^{4}) or higherorder accuracy is an artifact of the implementation in TempestRemap, which can be improved with further analysis and is not a limitation of the underlying numerical method.
TempestRemap produces conservative solution projections between mesh combinations that can be thirdorder with p=2 for smooth fields. However, even if local elementwise conservation can be guaranteed in overlaybased highorder 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 highorder linear maps that are also monotone. Due to this restriction, and since Gibbs phenomena (Jerri, 2013) are ubiquitous with highorder maps when steep gradients are encountered, global bounds can be preserved in TempestRemap by employing nonlinear filtering algorithms such as CAAS during the online solution transfer in ESMs.
5.1.3 GMLS
The remapping scheme based on the meshless generalized moving least squares (GMLS) method demonstrates the flexibility to deliver higherorder convergence for scalar field data. The convergence rates computed for the AnalyticalFun1 field are shown in Table 4 for various polynomial degrees.
The convergence rates for the nominal GMLS scheme and the GMLSCAAS remapping method with a postprocessing step in Table 4 show that, in general, highorder accuracy can be achieved. However, using the nominal GMLS scheme for climate modeling problems can result in nonconservative and potentially oscillatory reconstructions for fields with strong gradients. Hence, GMLSCAAS algorithm is especially advantageous to enable global and local bounds preservation. Note that the augmented GMLSCAAS algorithm suffers from convergence degradation for higher polynomial degree values and is limited to 𝒪(h^{3}) for this smooth field data in ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$. There is limited theoretical proof for convergence rates of the CAAS filter in the literature (Bradley et al., 2019) for arbitrary problems. However, in two dimensions, one can consider a field with 1D connected bands of extrema that demonstrates a maximum rate of convergence of 𝒪(h^{3}) in ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ norm. Algorithms like CAAS function by clipping newly formed local extrema resulting from higherorder reconstructions and computing a redistribution of the mass deficit accordingly. We hypothesize that the observed convergence degradation is primarily a result of these clipping and redistribution steps.
Even though the GMLSCAAS remaps show lower convergence rates, it provides the benefit of making the scheme globally conservative and monotone. The CAAS algorithm requires runtime modification of the projected fields to ensure global and local bounds, and the nonlinear solutiondependent filter can eliminate Gibbs oscillations, providing better stability during remap operation.
5.1.4 WLSENOR
The convergence rates for various polynomial degrees of reconstruction p are tabulated in Table 5. The convergence analysis for the WLSENOR scheme shows that even for high polynomial degrees, theoretically expected rates of 𝒪(h^{p+1}) in ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ and 𝒪(h^{p}) in ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ are observed for the smooth analytical fields. For example, for p=4, we can confirm that the asymptotic convergence rate for ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ is 𝒪(h^{5}) and for ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$, ${\leftE\right}_{{H}_{\mathrm{1}}}$ is 𝒪(h^{4}).
We note that the WLSENOR 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 highorder reconstruction process locally. In contrast to the GMLSCAAS highorder meshless scheme with a postprocessing filter that results in a convergence order degradation, the WLSENOR scheme remains consistently accurate for smooth functions up to p=4 in our experiments.
5.1.5 Real fields: convergence rate comparisons
While highorder convergence rates are achievable for smooth field profiles of AnalyticalFun1 and AnalyticalFun2, maintaining theoretical rates of convergence for “realworld” 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 𝒞^{0} and 𝒞^{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 SPHsampled TPW field data was computed, as tabulated in Table 6. The convergence rate comparisons using different remap schemes show that the computed rates for higherorder polynomial reconstructions fall severely short of theoretical rates as expected relative to the rates shown previously for AnalyticalFun1.
However, it can be observed that higher than firstorder 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 ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ and ${\leftE\right}_{{H}_{\mathrm{1}}}$. The meshless schemes are competitive in terms of accuracy bounds and provide a viable alternative approach for usage in production climate simulations, in which L^{2}minimization schemes and loworder bilinear maps have traditionally 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, p=1 or p=2 may typically yield a more numerically stable solution due to Gibbs phenomena. In such operational circumstances, highorder methods for TempestRemap and GMLS should be augmented with nonlinear postprocessing filters, such as CAAS.
5.2 Global conservation
All meshbased 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 sum operations on the discrete linear map matrix representing the projection or massmatrix operator.
The global field integral error indicating a conservation deficit for the TPW field variable is shown in Fig. 7 for different mesh resolutions of the CSMPAS mesh type. In these plots, the y axes represent the logarithmic error of the global integral value between the solution remapped onto Ω_{t} relative to the reference solution sampled on Ω_{t}. Hence, lower values indicate better satisfaction of the global conservation metric.
The expected behavior of preserving global integrals in the remapped solutions, as illustrated in Fig. 7, has been verified for both low and higherorder maps, irrespective of the mesh topology or the field variables being transferred across meshes. Notably, the accumulation of roundoff error in the global field integral as a function of remap iterations can quickly add to the global conservation error when using linear maps. This effect is seen in the field projections computed with both TempestRemap and ESMF through two sparse matrix–vector (SpMV) products in each remap iteration. As expected based on computational complexity of the SpMV operation, the loworder conservative ESMF remaps do accumulate errors more slowly than highorder TempestRemap projections, especially at fine mesh resolutions.
The WLSENOR scheme achieves excellent global conservation by applying adaptive quadrature rules to integrate the functions to (nearly) machine precision and redistributing the conservation errors locally near discontinuous regions. The unmodified GMLS scheme is nonconservative without any postprocessing and hence is not presented here. However, using the nonlinear CAAS filtering algorithm with the GMLS remapping scheme provides global conservation to userspecified tolerances.
This conservation metric describes whether the remapped solution preserves the global integral over the domain, which is often a desirable constraint for most scalar variables in order to ensure that the total mass and energy in the closed simulation system remain constant. However, for some scalar fields such as sea surface temperature (SST) or TPW, such a conservation constraint may not be strictly mandated, and this enables the use of even the nonconservative GMLS scheme, which has been demonstrated to achieve excellent accuracy.
5.3 Monotonicity – global extrema preservation
Departures away from global extrema, i.e., overshoots and undershoots, provide a useful metric to assess the stability of the remap operator under investigation. A perfectly stable method has a zero value for this metric. The results for the monotonicity metric as a function of the remap iterations are shown in Figs. 8 and 9 for $\left{G}_{\mathrm{max}}\right$ of the CFR variable on MPASRLL and $\left{G}_{\mathrm{min}}\right$ of the TOPO variable on CSMPAS meshes, respectively. These field variables have been specifically chosen to showcase the effects of the remapping algorithms on the welldefined upper and lower global field bounds that cannot be violated.
The results clearly demonstrate that the meshbased remapping schemes such as those in ESMF and TempestRemap implementations do not adhere to the global extrema (both maxima and minima) after the first linear remap operator application. The magnitude of this departure is resolutiondependent, but the behavior is consistent in all cases tested. The use of loworder, dissipative linear maps from ESMF shows a slow drift away from global maxima, which are nearly recovered after several remaps (around 700 in Fig. 8). However, the global minima increase continuously within the bounds of the reference solution. Hence, the conserve2nd option in ESMF damps the global maxima and amplifies the global minima. While not currently pursued in the paper, further analysis of the spectral properties of the loworder forward and reverse maps can shed light on the dampening properties of these ESMF linear maps.
In contrast, the highorder TempestRemap solutions show a drift away from the bounds in every case. The mildly damped increase in the metric as a function of remap operator application signals the presence of spurious highfrequency modes in the linear map, which is more prevalent as the degree of reconstruction increases. Obtaining highorder remapped solutions in addition to preserving monotonicity in these schemes will require postprocessing 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 GMLSCAAS and hybrid WLSENOR remap schemes maintain the global bounds for most test cases and fields tested in this study. Augmenting the nonmonotone GMLS scheme with the CAAS bounds preservation algorithm shows excellent, stable behavior for all cases. WLSENOR shows relatively high compliance with actual bounds because its nonoscillatory 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, WLSENOR would not apply the adaptive weights and the global extremes may oscillate randomly at a magnitude comparable to the local discretization error, as seen in Fig. 8b. Note that the hidden traces of $\left{G}_{\mathrm{min}}\right$ and $\left{G}_{\mathrm{max}}\right$ for GMLSCAAS and WLSENOR in Figs. 9 and 8, respectively, are at the expected levels of zero (monotone reconstructions). This is evident for the GMLSCAAS and WLSENOR methods in many of the experiments.
5.4 Locality – local extrema preservation
Similar to the global bounds metrics that measure solution monotonicity, the local bounds metrics compute the norm of error resulting from the comparison of the reference sampled field data against the projected field in a onering local neighborhood. This metric provides insight into the preservation of local field features on repeated remap operator applications. The computed locality metrics for maxima and minima for various mesh resolutions and realworld fields are shown in Figs. 10 and 11.
These results showcase the fact that the highorder TempestRemap implementations perform quite well with minimal feature loss for all fields on all meshes. In contrast, the loworder ESMF linear maps have high dissipation growth that is dependent on both the resolution and remap iterations. For highresolution meshes, the amplification of the dissipation at every iteration can even be undamped and grows linearly with the ESMF conserve2nd option. This result implies that as the number of remap iterations increase, more information on local features is lost during every transfer for these cases.
For the WLSENOR scheme, which uses a hybrid strategy to detect local features and adaptive quadrature rules to resolve a high discrepancy in mesh resolutions, its preservation of local bounds exhibits a strong dependence on the spatial resolution of Ω_{s} and Ω_{t}. However, unlike ESMF, the slope of metric indicating departure away from local bounds remains near zero as remap iterations are increased, signifying that the reconstructions remain stable on repeated applications without any further feature losses, probably due to the use of adaptive quadrature rules.
In the meshless remapping category, the GMLS schemes utilize the CAAS algorithm to clip and conserve fields in a local subset of the neighborhood for each evaluation point. For GMLSCAAS, 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 onering 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 GMLSCAAS in keeping with the spirit of its application as a purely meshless technique.
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 highorder polynomial reconstructions. These conclusions are similar to those observed in the current intercomparison experiments shown in Figs. 10 and 11. Note that in Fig. 10c, for equalresolution CSMPAS studies, the dissipation decreases with ESMF, GMLSCAAS, WLSENOR, and TempestRemap schemes in that order.
One of the key outcomes of this analysis indicates that using nonlinear filtering algorithms like CAAS can provide the benefits of property preservation to achieve global conservation and monotonicity constraints. However, these advantages can be offset by the higher local dissipation effects, especially on disparate mesh resolutions of Ω_{s} and Ω_{t}. So depending on the problem use case for which linear maps are being used to transfer solution fields between meshes with varying resolutions, an adaptive approach may be taken to choose a nonlinear filter when monotonicity constraints are important and to apply it with care, especially for highorder maps on which feature preservation may suffer.
5.5 Analysis on RRMs
We restrict the analysis in this study to all the conservative variations of the schemes chosen. As a result, the global conservation metric is truly satisfied for all cases tested with RRMs. The following subsections provide a detailed comparison of the error resolutions as well as global and local bound preservation metrics.
5.5.1 Error metric comparison
In order to compare the performance of the different remapping schemes on more realistic and complex unstructured RRMs, the global error norms ${\u2225E\u2225}_{{L}_{\mathrm{1}}}$, ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$, ${\u2225E\u2225}_{{L}_{\mathrm{\infty}}}$, ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$, and ${\leftE\right}_{{H}_{\mathrm{1}}}$ are computed and tabulated in Table 7. These results were computed on the CS(r2)–MPAS(r2) RRM combination to observe the accurate resolution of different remapped fields.
We note that the global errors with respect to all error norms are considerably smaller in WLSENOR and TempestRemap for the smooth (analytical) field variables sampled on the RRM input meshes. Interestingly, the loworder ESMF implementation shows performance comparable to the highorder GMLSCAAS(p=4) scheme in ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$, even though both these methods are more than 5 orders of magnitude worse than WLSENOR in terms of absolute error value. The gradient error metrics (${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ and ${\leftE\right}_{{H}_{\mathrm{1}}}$) also show the high dissipative errors in loworder ESMF methods for these smooth field data.
However, it is imperative to note that when transferring field data with significant C^{1} discontinuities (TOPO, CFR, TPW) that are typical in realworld climate simulations, using highorder remapping schemes such as WLSENOR, GMLS meshless methods and even TempestRemap(p=3) only provide nominal improvements over linear maps computed with loworder ESMF reconstructions with the conserve2nd option. These results demonstrate that the stateoftheart ESMF conservative remapping schemes can be complemented by other algorithms that have relatively better accuracy profiles and lower feature dissipation for many coupled climate modeling problems.
In order to make a fair selection of ESMF conservative algorithm, the comparison between the conserve and conserve2nd implementations was repeated on RRMs. The ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$ error profiles for various equalresolution source and target meshes are shown in Fig. 12 for the TPW variable. The results suggest that the behavior of the conserve2nd option remains marginally better than the conserve option in these experiments, even though the magnitude of accuracy gain with conserve2nd is not as attractive in comparison to the case of regularly refined meshes analyzed in Fig. 6. Hence, as an outcome of this experiment, the recommendation is to use conserve2nd in all existing production climate simulations when highorder maps are unavailable, as it has demonstrated more favorable properties than the conserve option in ESMF.
5.5.2 Monotonicity metrics
Similar to the uniformly refined cases, the monotonicity metric $\left{G}_{\mathrm{max}}\right$ was analyzed on the CS(r2)–MPAS(r2) regionally refined meshes for the CFR field, as plotted in Fig. 13. The behavior of the $\left{G}_{\mathrm{max}}\right$ metric shows similar trends in the RRM cases compared to the uniform refinement cases. The highorder TempestRemap method shows the worst behavior, while loworder ESMF maps exhibit good dampening after the first remap step. Note that the hidden traces of $\left{G}_{\mathrm{max}}\right$ for GMLSCAAS and WLSENOR in Fig. 13 are at the expected levels of zero (monotone reconstructions). The observation shown here confirms that overlaymeshbased schemes like ESMF and TempestRemap require additional postprocessing monotone filters such as CAAS during the online solution projection step in a simulation in order to ensure strict global bounds preservation and to avoid Gibbs oscillations for strong gradient fields.
5.5.3 Locality metrics
The observed trends in the locality metrics computed on RRMs are similar to those exhibited in the uniformly refined mesh experiments presented earlier. Hence, for brevity, only the ${L}_{max,\mathrm{2}}$ and L_{min,2} metrics are presented for the TOPO field on the finest resolution (r2) of RRMs in Fig. 14.
It is imperative to recognize that the highorder remapping schemes from TempestRemap and WLSENOR continue to generate minimally diffusive projections on target RRMs. This behavior is quite attractive as local features in the fields, even in the presence of strong gradients, are resolved accurately with very low dissipation away from the sampled reference field data. The locally adaptive, discontinuity tracking WLSENOR algorithm shows strongly stable behavior for all test cases.
However, the meshless GMLSCAAS scheme and the loworder ESMF scheme exhibit severe departures away from these local bounds that are consistent with the relatively larger ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ and ${\leftE\right}_{{H}_{\mathrm{1}}}$ observed in Table 7. Depending on the solution fields being remapped between model components, such strong dissipation in sharp features could lead to local numerical artifacts in the coupled climate system.
One key observation from the global monotonicity metrics shown in Fig. 13 is that TempestRemap fails to maintain strict global bounds, while loworder ESMF and the meshless GMLSCAAS schemes recover them consistently. However, the local bound preservation metrics in Fig. 14 demonstrate that when using the conservative, highorder TempestRemap algorithm, elementwise dissipation in the field data is strictly bounded and relatively much smaller in comparison to other remapping algorithms.
This work provides a foundation for the systematic computation of key metrics of interest in regridding problems and implements the infrastructure to sample analytical and realworld fields as well as to measure the metrics through a flexible Python code. In the current paper, we have presented detailed experiments to compare remapping schemes for FV–FV discretization settings on three specific mesh families (CS, MPAS, RLL) of varying resolutions, both uniform and regionally refined. There are several research directions that will be explored further in order to strengthen and improve the intercomparison protocol and apply it to more realistic climate science problems.
Four key extensions of interest in this context are

complex meshes, including successive radapted meshes, and realistic samples with topological holes;

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

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

computational efficacy, comparing numerical accuracy and algorithmic performance on nextgeneration architectures.
Further details regarding these research directions are provided in the following subsections.
6.1 Complex meshes: more than just global meshes
Uniformly refined meshes provide a simple infrastructure to test the consistency of remapping schemes, and additionally, regionally refined global grids help identify any spatialscale 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 highaltitude points from the Andes can influence and contaminate remaps to produce biases that appear along these coastal regions. Such regions can yield incorrect field remaps when using simple nearest neighbor maps, highorder reconstructions, or even when computing the overlay mesh, as the numerical tolerances for querying neighborhood data become an important influencing factor in the overall accuracy.
The metrics definitions for standard error norms, global maxima, and global minima introduced in Sect. 3 are valid for general meshes. But, in order for the global conservation metric L_{g} shown in Eq. (11) to be valid, the target domain has to fully envelop the source domain. Additionally, the definition and implementation of local dissipation metrics use a small neighborhood (onering patch) around the target cell, which would also impose a similar requirement. This implies that for analysis with realistic meshes, the metrics for evaluating the quality of solution field data transferred between the atmosphere and ocean components will strictly be measured only on the global atmosphere meshes to satisfy the above constraints. Note that the input specification for such cases would require at least the area fractions and appropriate masks to be provided in order to enable efficient remap computations using the meshless and hybrid schemes like GMLS and WLSENOR.
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 “cyclicremap” test suites in the ALE community (Bochev and Shashkov, 2005). Performing successive remap over the sequence of such meshes, one can directly compare the final resulting field after multiple remaps directly with the initial field without a need for a reference analytical solution field. Devising such a study requires 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
The current study focused primarily on elementaveraged field data typical of FV discretization of models. Although many atmosphere and ocean models define the scalar fields as elementaveraged 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 highorder, conservative remapping schemes tested in the current study under discretization specifications commonly used in climate components like HOMME (Taylor et al., 2007). We also make a note that such extensions to more complex discretizations of source and target data should not require any modifications in the metrics definitions already presented in this paper. Further extensions to meshes that contain topological holes will still require Ω_{s}⊂Ω_{t}.
6.3 Complex fields: more than just scalars
During the remapping of climate data in production runs and scientific analysis, it is important to preserve not only scalar fields but also vector fields and derived properties that provide better insights from simulations. For example, to understand the atmospheric flows, analysis of global wind patterns is often performed in weather prediction, which requires the evaluation of derivatives or integrals of the vector wind velocity fields (Nair and Jablonowski, 2008). Some critical considerations in such scenarios require preservation of divergencefree velocity fields, conservation of vorticity, and preservation of windstress 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 divergencefree conditions and can be inconsistent, since the propagation of remapping error in the components is not correlated. Care should also be taken to ensure that the regridding of vector components uses a proper 3D Cartesian coordinate system instead of spherical mesh projections, which will not yield consistent vector fields on Ω_{t}.
Remapping algorithms based on mimetic schemes (Pletzer and Hayek, 2019) that provide exact conservation for both scalar and vector fields are promising in this direction. To our knowledge, existing remapping algorithms based on L^{2} minimization and highorder interpolationbased algorithms need further research for tackling vector field data in problems of interest. Furthermore, field data with nonlinear correlations such as those analyzed by Lauritzen and Thuburn (2012), adapted for remapping scenarios, can also be valuable in determining whether the solution transfer algorithms can remain conservative and preserve correlation properties. Nonlinear remapping schemes (Carey et al., 2001; Bochev and Shashkov, 2005; Bochev et al., 2011) may also prove to be viable options for these cases, although computational cost can be relatively much higher than using linear maps on decoupled scalar components. Additionally, some remapping metrics definitions introduced in Sect. 3 will have to be extended for vector field data.
6.4 Computational efficacy: comparing accuracy and efficiency
The additional characteristic dimension that is essential to recognize in intercomparison studies is that the overall cost to obtain the remaps over the simulation cycle includes both a onetime setup cost and a constantly growing computational load at every coupling step when field data need to be transferred between components. Hence, in order to better describe the computational cost of the remapping algorithms, we could split the effort into an offline cost and an online cost. Note that, often, what is sacrificed in terms of computational performance is gained in the quality of the remapped solution through higher accuracy, global bounds preservation, and strong local feature resolution with minimal dissipation. Since the numerical efficiency and time to compute the solution are competing factors, the usage of advanced remapping algorithms for realistic cases will require a more detailed analysis of the computational complexity at scale. An exception to the performance model derived here is that for dynamic or moving meshes such as those encountered in sea ice and ocean coupling in ESMs, the offline cost is nearly zero and the entire computational cost is online, since the maps to compute field projections have to be recomputed with changes in Ω_{s} and Ω_{t}.
The traditional meshbased conservative remapping algorithms used in ESMF and TempestRemap require computation of an intersection or supermesh, Ω_{s∩t}, which in general is computationally expensive. While linear complexity strategies like advancing front methods do exist in libraries like PANG (Gander and Japhet, 2013) and MOAB (Mahadevan et al., 2020) to compute mesh intersections, they can suffer from robustness issues when Ω_{s} has topological holes. An alternative approach is to utilize a K–d tree data structure to accelerate queries on the meshes in order to exactly locate the element containing a particular point, which is a fundamental operation in intersection mesh computation (Jansen et al., 1992). The leading computational cost for this operation is 𝒪(Nlog (N)), where N is the total number of elements in the mesh. The linear map computation itself follows typical FV or FE operator assembly workflows and remains bounded. The online cost of such a linear map application is essentially then provided by the sparse matrix–vector (SpMV) operation. The complexity for SpMV is, however, dependent on the degree p used for field reconstruction as it dictates the bandwidth of coupling between degrees of freedom in Ω_{s} and Ω_{t}. The illustration in Fig. 15 for ESMF (conserve) and TR (p=3) remapping algorithms shows the number of CS elements required to reconstruct the solution on a single target MPAS element when computing projections between the CS(r3)–MPAS(r3) uniform refinement case. While the loworder ESMF algorithm may require only 2 FLOPS per reconstruction on a target element, TempestRemap (p=3) would require a minimum of 50 FLOPS to achieve better accuracy. To overcome some of these computational issues for large meshes, the use of MPI parallelism has been exploited in ESMF and has been proven to be scalable with MOAB–TempestRemap (Mahadevan et al., 2020) libraries.
Computation for the GMLSCAAS meshfree approach on a manifold is dominated by the internal QR with pivoting factorization, requiring $\mathcal{O}\left(\frac{({p}^{\mathrm{2}}{)}^{\mathrm{3}}}{\mathrm{6}}\cdot N\right)$ FLOPS for the offline stencil calculation, whereby p is the polynomial degree of the basis used. This stencil calculation can be stored and applied as an SpMV operation, similar to TempestRemap. The CAAS nonlinear filter is an online cost as it is solutiondependent, so at every instance that it is applied it carries a $\mathcal{O}\left(\frac{{p}^{\mathrm{2}}}{\mathrm{2}}\right)$ bounds calculation cost for each of N degrees of freedom (DoF).
In terms of the computational cost, similar to TempestRemap and GMLS, the transfer for smooth function and the detection of discontinuities in WLSENOR can take advantage of some preprocessing steps to build matrixbased operators. Hence, they primarily involve SpMV as the primary online cost as well. Additionally, the resolution of discontinuities requires constructing and solving the generalized Vandermonde systems for each target cell near discontinuities as described in Sect. 2.2, which can be relatively expensive. Hence, the overall cost for resolving discontinuities depends on the percentage of the cells in discontinuity regions, which is determined only at runtime to account for field evolution. In general, on coarse meshes, this ratio may be relatively higher. As the underlying mesh is refined, the solution field profile has a finer approximation that results in a lower proportion of discontinuous cells. So the overall cost for the most expensive part in WLSENOR is expected to have a sublinear time complexity under mesh refinement.
While the theoretical complexity requirements for the remapping algorithms under consideration can be considerably different, it is imperative to measure the performance of these implementations on various architectures on standard test problems to gauge overall efficacy (accuracy vs. total computational time for N_{r} steps). This requires configuring, building, and installing the ESMF, TempestRemap, Compadre (GMLSCAAS), and WLSENOR libraries on the same machine to compare the performance profile for test problems at scale. Given the complexity and magnitude of such a task, we express this as another future avenue for research experiments that can add value to the broader climate science community using remapping algorithms (Valcke et al., 2022).
Remapping algorithms are critically important in climate science applications to maintain key numerical properties during the transformation and transfer of coupled field data between component models. Inaccurate, nonconservative, or highly dissipative remapping operators can introduce numerical artifacts in the coupled fields, which can destroy the highorder accuracy of component solvers and propagate errors into the global nonlinear system representing the climate system. Hence, understanding the behavior of remapping algorithms for Earth system modeling requires standardized numerical definitions that provide better insight into the accuracy of transferred fields, along with the ability to preserve global and local solution bounds to avoid numerically induced instabilities in the system.
In this paper, with these motivations in mind, four different remapping schemes were selected based on the current availability of the software implementations, the maturity of the underlying numerics, and potential computational efficiency that could be obtained for realworld scenarios in comparison to existing stateoftheart remapping algorithms used in production climate simulations. The meshbased implementations in the libraries ESMF (Hill et al., 2004) and TempestRemap (Ullrich and Taylor, 2015; Ullrich et al., 2016) provide support for the projection of field data between unstructured grids by using intersection meshes to compute conservative weights through global L^{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 intersectionmeshbased remapping approaches, a hybrid meshbased approach without an explicit need for overlay meshes using the WLSENOR 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 postprocessing filters to ensure monotonicity in the projected solutions using the ClipAndAssuredSum (CAAS) (Bradley et al., 2019) technique. These remapping algorithms cover a large span of low and highorder conservative solution transfer implementations that can directly impact the overall stability and accuracy of predictive solvers and analysis suites for weather and climate modeling.
Comparing these four distinct remapping algorithms requires a uniform test infrastructure, which provides the framework to create new verification studies and analyze the metrics obtained from the remapping algorithms. To enable such unbiased comparative studies, we introduced several remapping metrics that represent the key properties of remapping algorithms. These include global error measures under various norms like ${\u2225E\u2225}_{{L}_{\mathrm{1}}}$, ${\u2225E\u2225}_{{L}_{\mathrm{2}}}$, and ${\u2225E\u2225}_{{L}_{\mathrm{\infty}}}$ and gradient error measures given by ${\u2225E\u2225}_{{H}_{\mathrm{1}}}$ and ${\leftE\right}_{{H}_{\mathrm{1}}}$. These error norms provide the necessary verification of the theory and implementation of the remapping algorithms by measuring the theoretical order of accuracy applied to solution fields with sufficient smoothness. Next, global conservation errors are measured by comparing the integral of the sampled reference fields on Ω_{t} against the remapped solution on Ω_{t} for multiple projection iterations. Finally, metrics for the global and local departure away from maxima or minima in the solution fields can be evaluated in various norms to provide insight on monotonicity preservation and feature dissipation due to the remapping algorithms. These standardized metrics provide the blueprint to build the remapping intercomparison suite, which was then used to understand the numerical properties of all algorithms under consideration.
Furthermore, a flexible workflow built on several Pythonbased drivers to generate the unstructured meshes of different element topologies and resolutions, including regionally refined meshes, and to accurately sample five elementaveraged fields using SPH expansions was provided. With these input meshes, the four remapping algorithms were applied for both the smooth analytical fields and representative real fields in an iterative fashion to compute cyclic projections (${\mathrm{\Omega}}_{\mathrm{s}}\to {\mathrm{\Omega}}_{\mathrm{t}}\to {\mathrm{\Omega}}_{\mathrm{s}}$) for FV–FV field transfers between component models in a climate system.
(Guerra et al., 2021)(Guerra et al., 2021)Guerra and Mahadevan (2021)(Wimmers and Velden, 2011)(Platnick et al., 2020)(NOAA National Geophysical Data Center, 2009; Amante and Eakins, 2009)(Guerra and Mahadevan, 2021)Mahadevan et al. (2021)(Mahadevan et al., 2021)The results compiled from various test problems demonstrate that the conservative remapping implementations in ESMF are all firstorder accurate for smooth problems, even though the conserve2nd option produces consistently better accuracy than the conserve option. Hence, when possible, the ESMF conserve2nd option should be used to obtain better remaps for climate simulations. However, these loworder ESMF maps are highly dissipative in general. In contrast, TempestRemap produces remaps that are globally conservative with highorder accurate convergence to 𝒪(h^{p+1}) for smooth solution fields using seconddegree (p=2) polynomial reconstructions and bounded dissipation on field projections containing sharp features. Here, h represents the characteristic spatial length of the mesh wherein errors are measured. However, neither of these low and highorder L^{2}minimization approaches can guarantee monotone remaps without the use of postprocessing filters such as slope limiters or CAAS to enforce global and local solution bounds. The use of these filters, however, inevitably introduces additional dissipation that can become significant with repeated applications of the remap operator.
On the other hand, the hybrid meshbased and meshless schemes achieve very high order consistently (up to 𝒪(h^{5})) while retaining global conservation, which has been one of the key advantages of traditional overlapmeshbased schemes implemented in ESMF and TempestRemap. Additionally, the ability to capture smoothly varying fields very accurately without any degradation even for highpolynomialdegree reconstructions makes these schemes attractive and competitive compared to traditional meshbased 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 be especially valuable when highly accurate scalar fields need to be sampled on a refined RLL grid for further analysis or in situ visualization to track multidecadal climate evolution.
It is essential to note that while the fully meshbased remapping algorithms are, in general, insensitive to mesh resolutions or the topology, both the hybrid and meshless schemes are susceptible to larger meshdependent dissipation. This is especially evident when the source and target mesh resolutions differ drastically. The use of the CAAS algorithm, when combined with even nonconservative schemes such as GMLS, can provide global conservation and bounds preservation at the cost of added dispersion to control numerical oscillations. In contrast, the builtin discontinuity indicators used by the WLSENOR algorithm demonstrate good featureresolving properties for real fields in all experiments conducted.
These experiments conducted on both uniformresolution 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 realworld scenarios requires deeper investigation using more topologically diverse, complex meshes and discretization specifications that are more representative of components used in E3SM and CESM.
Finally, we want to emphasize that the MIRA infrastructure presented in this paper is freely available as an opensource package (Guerra et al., 2021) to compare new and existing remapping algorithms under the same overall test constraints. Such intercomparison studies are important to evaluate the cost of remapping algorithms under stability and accuracy constraints, which remain crucial to better understanding the propagation of errors in coupled climate and weather systems.
Information on the availability of source code for the remapping metrics intercomparison infrastructure featured in this paper, all relevant input meshes, and the final consolidated metrics data for schemes are provided in Table 8.
VSM, JEG, XJ, and PK wrote the paper (with several helpful review comments from PB, PU, and RJ). DM contributed text to the TempestRemap section. JG wrote significant portions of the remapping intercomparison code (Guerra et al., 2021), with contributions from VSM and PK. VSM generated the remap output data using ESMF and TempestRemap libraries. PK generated results using GMLS and GMLSCAAS. XJ and YL generated results using WLSENOR algorithms. VSM 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.
At least one of the (co)authors is a member of the editorial board of Geoscientific Model Development. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We gratefully acknowledge the computing resources provided on Bebop, a highperformance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract DEAC0206CH11357, 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 WLSENOR were obtained using the Seawulf cluster at the Institute for Advanced Computational Science of Stony Brook University, which was partially funded by Empire State Development grant NYS no. 28451. The authors would also like to thank the two anonymous reviewers for their detailed comments, which improved the clarity of a few sections in this paper.
This research is supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (grant no. KJ0403000) as well as the Office of Biological and Environmental Research (grant no. KP1703020) Scientific Discovery through Advanced Computing (SciDAC) program. Additional support was provided by the NOAA Office of Oceanic and Atmospheric Research under the NOAA–University of Oklahoma cooperative agreement (NA11OAR4320072; U.S. Department of Commerce).
This paper was edited by Simone Marras and reviewed by two anonymous referees.
Amante, C. and Eakins, B. W.: ETOPO1 1 ArcMinute Global Relief Model: Procedures, Data Sources and Analysis, NOAA Technical Memorandum NESDIS NGDC24, https://doi.org/10.7289/V5C8276M, 2009. a, b
Balaji, V., Anderson, J., Held, I., Winton, M., Durachta, J., Malyshev, S., and Stouffer, R. J.: The Exchange Grid: A mechanism for data exchange between Earth System components on independent grids, in: Parallel Computational Fluid Dynamics 2005, edited by: Deane, A., Ecer, A., McDonough, J., Satofuka, N., Brenner, G., Emerson, D. R., Periaux, J., and TromeurDervout, D., Elsevier, Amsterdam, 179–186, https://doi.org/10.1016/B9780444522061/500215, 2006. a
Barth, T. and Jespersen, D.: The design and application of upwind schemes on unstructured meshes, in: 27th Aerospace sciences meeting, 366 pp., https://doi.org/10.2514/6.1989366, 1989. a
Berger, M., Murman, S. M., and Aftosmis, M. J.: Analysis of Slope Limiters on Irregular Grids, in: Proceedings of the 43rd AIAA Aerospace Sciences Meeting, AIAA20050490, AIAA, Reno, NV, 2005. a
Blanchard, G. and Loubere, R.: HighOrder Conservative Remapping with a posteriori MOOD stabilization on polygonal meshes, Comput. Fluids, 136, 83–103, https://doi.org/10.1016/j.compfluid.2016.06.002, 2016. a
Bochev, P. and Shashkov, M.: Constrained interpolation (remap) of divergencefree fields, Comput. Meth. Appl. Mech., 194, 511–530, 2005. a, b, c
Bochev, P., Ridzal, D., Scovazzi, G., and Shashkov, M.: Formulation, analysis and numerical study of an optimizationbased conservative interpolation (remap) of scalar fields for arbitrary Lagrangian–Eulerian methods, J. Comput. Phys., 230, 5199–5225, 2011. a, b
Bochev, P., Ridzal, D., and Peterson, K.: Optimizationbased remap and transport: A divide and conquer strategy for featurepreserving discretizations, J. Comput. Phys., 257, 1113–1139, 2014. a, b
Bradley, A. M., Bosler, P. A., Guba, O., Taylor, M. A., and Barnett, G. A.: Communicationefficient property preservation in tracer transport, SIAM J. Sci. Comput., 41, C161–C193, 2019. a, b, c
Breitkopf, P., Rassineux, A., Touzot, G., and Villon, P.: Explicit form and efficient computation of MLS shape functions and their derivatives, Int. J. Numer. Meth. Eng., 48, 451–466, 2000. a
Brewer, M. L., Diachin, L. F., Knupp, P. M., Leurent, T., and Melander, D. J.: The Mesquite Mesh Quality Improvement Toolkit, in: Proceedings of the 12th International Meshing Roundtable, IMR 2003, Santa Fe, New Mexico, USA, 14–17 September 2003. a
Buhmann, M.: A new class of radial basis functions with compact support, Math. Comput., 70, 307–318, 2001. a
Bungartz, H.J., Lindner, F., Gatzhammer, B., Mehl, M., Scheufele, K., Shukaev, A., and Uekermann, B.: preCICE – A fully parallel library for multiphysics surface coupling, Comput. Fluids, 141, 250–258, 2016. a, b, c
Carey, G., Bicken, G., Carey, V., Berger, C., and Sanchez, J.: Locally constrained projections on grids, Int. J. Numer. Meth. Eng., 50, 549–577, 2001. a, b
NOAA National Geophysical Data Center: ETOPO1 1 ArcMinute Global Relief Model, NOAA National Centers for Environmental Information, 2009. a, b
Chesshire, G. and Henshaw, W. D.: A scheme for conservative interpolation on overlapping grids, SIAM J. Sci. Comput., 15, 819–845, 1994. a
Collins, N., Theurich, G., DeLuca, C., Suarez, M., Trayanov, A., Balaji, V., Li, P., Yang, W., Hill, C., and Da Silva, A.: Design and implementation of components in the Earth System Modeling Framework, Int. J. High Perform. C., 19, 341–350, 2005. a
Craig, A., Valcke, S., and Coquart, L.: Development and performance of a new version of the OASIS coupler, OASIS3MCT_3.0, Geosci. Model Dev., 10, 3297–3308, https://doi.org/10.5194/gmd1032972017, 2017. a, b
Craig, A. P., Vertenstein, M., and Jacob, R.: A new flexible coupler for earth system modeling developed for CCSM4 and CESM1, Int. J. High Perform. C., 26, 31–42, https://doi.org/10.1177/1094342011428141, 2012. a
de Boer, A., van Zuijlen, A. H., and Bijl, H.: Comparison of conservative and consistent approaches for the coupling of nonmatching meshes, Comput. Method. Appl. M., 197, 4284–4297, 2008. a, b
de Boor, C.: Quasiinterpolants and approximation power of multivariate splines, in: Computation of curves and surfaces, Springer, 313–345, 1990. a
Dukowicz, J. K. and Baumgardner, J. R.: Incremental Remapping as a Transport/Advection Algorithm, J. Comput. Phys., 160, 318– 335, https://doi.org/10.1006/jcph.2000.6465, 2000. a
Dukowicz, J. K. and Kodis, J. W.: Accurate conservative remapping (rezoning) for arbitrary LagrangianEulerian computations, SIAM J. Sci. Stat. Comp., 8, 305–321, 1987. a, b, c, d
E3SM Project: Energy Exascale Earth System Model (E3SM), https://doi.org/10.11578/E3SM/dc.20180418.36, 2018. a, b
Edwards, H. C., Trott, C. R., and Sunderland, D.: Kokkos: Enabling manycore performance portability through polymorphic memory access patterns, J. Parallel Distr. Com., 74, 3202–3216, https://doi.org/10.1016/j.jpdc.2014.07.003, 2014. a
Erath, C., Lauritzen, P. H., and Tufo, H. M.: On mass conservation in highorder highresolution rigorous remapping schemes on the sphere, Mon Weather Rev., 141, 2128–2133, 2013. a, b
Farrell, P. and Maddison, J.: Conservative interpolation between volume meshes by local Galerkin projection, Comput. Method. Appl. M., 200, 89–100, 2011. a
Farrell, P., Piggott, M., Pain, C., Gorman, G., and Wilson, C.: Conservative interpolation between unstructured meshes via supermesh construction, Comput. Method. Appl. M., 198, 2632–2642, 2009. a, b, c
Flyer, N. and Wright, G. B.: Transport schemes on a sphere using radial basis functions, J. Comput. Phys., 226, 1059–1084, 2007. a
Fornberg, B. and Flyer, N.: The Gibbs phenomenon for radial basis functions, in: The Gibbs Phenomenon in Various Representations and Applications, Potsdam, NY: Sampling Publishing, 201–224, 2007. a
Gander, M. J. and Japhet, C.: Algorithm 932: PANG: software for nonmatching grid projections in 2D and 3D with linear complexity, ACM T. Math. Softw., 40, 1–25, 2013. a, b, c
Garimella, R., Kucharik, M., and Shashkov, M.: An efficient linearity and bound preserving conservative interpolation (remapping) on polyhedral meshes, Comput. Fluids, 36, 224–237, 2007. a
Godunov, S. K.: A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Mat. Sb., 89, 271–306, 1959. a, b
Golub, G. H. and Van Loan, C. F.: Matrix Computations, Johns Hopkins, 4th Edn., ISBN 1421407949 9781421407944, 2013. a, b
Gottlieb, D. and Shu, C.W.: On the Gibbs phenomenon and its resolution, SIAM Rev., 39, 644–668, 1997. a
Grandy, J.: Conservative remapping and region overlays by intersecting arbitrary polyhedra, J. Comput. Phys., 148, 433–466, 1999. a, b
Gross, B., Trask, N., Kuberry, P., and Atzberger, P.: Meshfree methods on manifolds for hydrodynamic flows on curved surfaces: A Generalized Moving LeastSquares (GMLS) approach, J. Comput. Phys., 409, 109340, https://doi.org/10.1016/j.jcp.2020.109340, 2020. a, b, c
Guerra, J. and Mahadevan, V.: Satellite datasets used for MIRA workflows, Zenodo [data set], https://doi.org/10.5281/zenodo.5172792, 2021. a, b, c
Guerra, J., Mahadevan, V., Kuberry, P., Jiao, X., and Li, Y.: MIRA: Metrics for Intercomparison of Remapping Algorithms, Zenodo [code], https://doi.org/10.5281/zenodo.5518037, 2021. a, b, c, d, e
Hanke, M., Redler, R., Holfeld, T., and Yastremsky, M.: YAC 1.2.0: new aspects for coupling software in Earth system modelling, Geosci. Model Dev., 9, 2755–2769, https://doi.org/10.5194/gmd927552016, 2016. a, b, c, d
Hill, C., DeLuca, C., Balaji, Suarez, M., and Da Silva, A.: The architecture of the earth system modeling framework, Comput. Sci. Eng., 6, 18–28, https://doi.org/10.1109/MCISE.2004.1255817, 2004. a, b, c, d, e, f
Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., Lamarque, J.F., Large, W. G., Lawrence, D., Lindsay, K., Lipscomb, W. H., Long, M. C., Mahowald, N., Marsh, D. R., Neale, R. B., Rasch, P., Vavrus, S., Vertenstein, M., Bader, D., Collins, W. D., Hack, J. J., Kiehl, J., and Marshall, S.: The community earth system model: a framework for collaborative research, B. Am. Meterol. Soc., 94, 1339–1360, https://doi.org/10.1175/BAMSD1200121.1, 2013. a, b
Jalali, A. and Gooch, C. F. O.: HigherOrder Finite Volume Solution Reconstruction on Highly Anisotropic Meshes, in: 21st AIAA Computational Fluid Dynamics Conference, San Diego, https://doi.org/10.2514/6.20132565, 2013. a
Jansen, K., Shakib, F., and Hughes, T. J.: Fast projection algorithm for unstructured meshes, Computational nonlinear mechanics in aerospace engineering, 146, 175, 1992. a, b
Jerri, A. J.: The Gibbs Phenomenon in Fourier Analysis, Splines and Wavelet Approximations, Mathematics and Its Applications, Springer, 446, 340 pp., https://doi.org/10.1007/9781475728477, 2013. a, b
Jiao, X. and Heath, M. T.: Commonrefinementbased data transfer between nonmatching meshes in multiphysics simulations, Int. J. Numer. Meth. Eng., 61, 2402–2427, 2004a. a, b, c, d, e
Jiao, X. and Heath, M. T.: Overlaying surface meshes, part I: Algorithms, Int. J. Comput. Geom. Ap., 14, 379–402, 2004b. a
Jiao, X. and Heath, M. T.: Overlaying surface meshes, part II: Topology preservation and feature matching, Int. J. Comput. Geom. Ap., 14, 403–419, 2004c. a
Jiao, X. and Wang, D.: Reconstructing highorder surfaces for meshing, Eng. Comput., 28, 361–373, 2012. a
Joldes, G. R., Chowdhury, H. A., Wittek, A., Doyle, B., and Miller, K.: Modified moving least squares with polynomial bases for scattered data approximation, Appl. Math. Comput., 266, 893–902, 2015. a
Jones, P. W.: First and secondorder conservative remapping schemes for grids in spherical coordinates, Mon. Weather Rev., 127, 2204–2210, 1999. a, b, c, d, e, f, g, h, i, j, k
Joppich, W. and Kürschner, M.: MpCCI – a tool for the simulation of coupled applications, Concurr. Comp.Pract. E., 18, 183–192, 2006. a
Kageyama, A. and Sato, T.: “YinYang grid”: An overset grid in spherical geometry, Geochem. Geophy. Geosy., 5, Q09005, https://doi.org/10.1029/2004GC000734, 2004. a
Kritsikis, E., Aechtner, M., Meurdesoif, Y., and Dubos, T.: Conservative interpolation between general spherical meshes, Geosci. Model Dev., 10, 425–431, https://doi.org/10.5194/gmd104252017, 2017. a, b, c
Kuberry, P., Bosler, P., and Trask, N.: Compadre Toolkit, https://doi.org/10.11578/dc.20190411.1, 2019. a
Lancaster, P. and Salkauskas, K.: Surfaces generated by moving least squares methods, Math. Comput., 37, 141–158, 1981. a, b
Larson, J., Jacob, R., and Ong, E.: The model coupling toolkit: a new Fortran90 toolkit for building multiphysics parallel coupled models, Int. J. High Perform. C., 19, 277–292, 2005. a
Lauritzen, P. H. and Nair, R. D.: Monotone and conservative cascade remapping between spherical grids (CaRS): Regular latitude–longitude and cubedsphere grids, Mon. Weather Rev., 136, 1416–1432, 2008. a, b, c, d
Lauritzen, P. H. and Thuburn, J.: Evaluating advection/transport schemes using interrelated tracers, scatter plots and numerical mixing diagnostics, Q. J. Roy. Meteor. Soc., 138, 906–918, 2012. a
Lauritzen, P. H., Nair, R. D., and Ullrich, P. A.: A conservative semiLagrangian multitracer transport scheme (CSLAM) on the cubedsphere grid, J. Comput. Phys., 229, 1401–1424, 2010. a
Li, Y., Chen, Q., Wang, X., and Jiao, X.: WLSENO Remap: Superconvergent and NonOscillatory Weighted Least Squares Data Transfer on Surfaces, J. Comput. Phys., 417, 109578, https://doi.org/10.1016/j.jcp.2020.109578, 2020. a, b, c, d, e, f, g, h, i, j, k, l
Liang, J. and Zhao, H.: Solving partial differential equations on point clouds, SIAM J. Sci. Comput., 35, A1461–A1486, 2013. a, b
Liu, H. and Jiao, X.: WLSENO: Weightedleastsquares based essentially nonoscillatory schemes for finite volume methods on unstructured meshes, J. Comput. Phys., 314, 749–773, 2016. a, b
Liu, L., Zhang, C., Li, R., Wang, B., and Yang, G.: CCoupler2: a flexible and userfriendly community coupler for model coupling and nesting, Geosci. Model Dev., 11, 3557–3586, https://doi.org/10.5194/gmd1135572018, 2018. a, b
Mahadevan, V., Guerra, J., Kuberry, P., and Jiao, X.: MIRADatasets: Datasets from Metrics for Intercomparison of Remapping Algorithms, Zenodo [data set], https://doi.org/10.5281/zenodo.5518065, 2021. a, b, c, d, e
Mahadevan, V. S., Grindeanu, I., Jacob, R., and Sarich, J.: Improving climate model coupling through a complete mesh representation: a case study with E3SM (v1) and MOAB (v5.x), Geosci. Model Dev., 13, 2355–2377, https://doi.org/10.5194/gmd1323552020, 2020. a, b, c, d, e
Mirzaei, D., Schaback, R., and Dehghan, M.: On generalized moving least squares and diffuse derivatives, IMA J. Numer. Anal., 32, 983–1000, https://doi.org/10.1093/imanum/drr030, 2012. a
Nair, R. D. and Jablonowski, C.: Moving vortices on the sphere: A test case for horizontal advection problems, Mon. Weather Rev., 136, 699–711, 2008. a
Nayroles, B., Touzot, G., and Villon, P.: Generalizing the finite element method: Diffuse approximation and diffuse elements, Comput. Mech., 10, 307–318, https://doi.org/10.1007/BF00364252, 1992. a
Norman, M. R. and Nair, R. D.: Inherently conservative nonpolynomialbased remapping schemes: Application to semiLagrangian transport, Mon. Weather Rev., 136, 5044–5061, 2008. a
Petersen, M.: MPASOcean V6 Run Directories, Zenodo [code], https://doi.org/10.5281/zenodo.1252437, 2018. a
Platnick, S., Ackerman, S. A., King, M. D., Meyer, K., Menzel, W. P., Holz, R. E., Baum, B. A., and Yang, P.: MODIS atmosphere L2 cloud product (06_L2), NASA MODIS Adaptive Processing System [data set], https://doi.org/10.5067/MODIS/MYD06_L2.006, 2020. a, b
Pletzer, A. and Hayek, W.: Mimetic interpolation of vector fields on Arakawa C/D grids, Mon. Weather Rev., 147, 3–16, 2019. a, b
Rider, W. J.: Reconsidering remap methods, Int. J. Numer. Meth. Fl., 76, 587–610, https://doi.org/10.1002/fld.3950, 2014. a
Ringler, T. D., Thuburn, J., Klemp, J. B., and Skamarock, W. C.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarilystructured Cgrids, J. Comput. Phys., 229, 3065–3090, 2010. a
Royer, J.F.: Correction of negative mixing ratios in spectral models by global horizontal borrowing, Mon. Weather Rev., 114, 1406–1410, 1986. a
Skamarock, W. C. and Gassmann, A.: Conservative Transport Schemes for Spherical Geodesic Grids: HighOrder Flux Operators for ODEBased Time Integration, Mon. Weather Rev., 139, 2962–2975, https://doi.org/10.1175/MWRD1005056.1, 2011. a
Skamarock, W. C. and Menchaca, M.: Conservative Transport Schemes for Spherical Geodesic Grids: HighOrder Reconstructions for ForwardinTime Schemes, Mon. Weather Rev., 138, 4497–4508, https://doi.org/10.1175/2010MWR3390.1, 2010. a
Slattery, S., Wilson, P., and Pawlowski, R.: The Data Transfer Kit: A geometric rendezvousbased tool for multiphysics data transfer, in: International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering (M&C 2013), 5–9, ISBN 9780894487002, 2013. a, b
Slattery, S. R.: Meshfree data transfer algorithms for partitioned multiphysics problems: Conservation, accuracy, and parallelism, J. Comput. Phys., 307, 164–188, 2016. a, b, c, d
Smith, M. J., Cesnik, C. E., and Hodges, D. H.: Evaluation of some data transfer algorithms for noncontiguous meshes, J. Aerospace Eng., 13, 52–58, 2000. a
Suchde, P. and Kuhnert, J.: A meshfree generalized finite difference method for surface PDEs, Comput. Math. Appl., 78, 2789–2805, https://doi.org/10.1016/j.camwa.2019.04.030, 2019. a, b
Tautges, T. J. and Caceres, A.: Scalable parallel solution coupling for multiphysics reactor simulation, in: Journal of Physics: Conference Series, vol. 180, 012017, IOP Publishing, 2009. a
Taylor, M., Edwards, J., Thomas, S., and Nair, R.: A mass and energy conserving spectral element atmospheric dynamical core on the cubedsphere grid, J. Phys. Conf. Ser., 78, 012074, https://doi.org/10.1088/17426596/78/1/012074, 2007. a
Thuburn, J., Ringler, T. D., Skamarock, W. C., and Klemp, J. B.: Numerical representation of geostrophic modes on arbitrarily structured Cgrids, J. Comput. Phys., 228, 8321–8335, 2009. a
Townsend, A., Wilber, H., and Wright, G. B.: Computing with functions in spherical and polar geometries I. The sphere, SIAM J. Sci. Comput., 38, C403–C425, 2016. a
Trask, N. and Kuberry, P.: Compatible meshfree discretization of surface PDEs, Comput. Particle Mech., 7, 271–277, https://doi.org/10.1007/s40571019002512, 2020. a, b, c, d, e
Ullrich, P. A. and Taylor, M. A.: Arbitraryorder conservative and consistent remapping and a theory of linear maps: Part I, Mon. Weather Rev., 143, 2419–2440, 2015. a, b, c, d, e, f, g, h, i, j
Ullrich, P. A., Lauritzen, P. H., and Jablonowski, C.: Geometrically Exact Conservative Remapping (GECoRe): Regular latitude–longitude and cubedsphere grids, Mon. Weather Rev., 137, 1721–1741, 2009. a, b, c
Ullrich, P. A., Devendran, D., and Johansen, H.: Arbitraryorder conservative and consistent remapping and a theory of linear maps: Part II, Mon. Weather Rev., 144, 1529–1549, 2016. a, b, c, d
Valcke, S., Piacentini, A., and Jonville, G.: Benchmarking Regridding Libraries Used in Earth System Modelling, Math. Comput. Appl., 27, https://doi.org/10.3390/mca27020031, 2022. a, b
Van Leer, B.: Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov's method, J. Comput. Phys., 32, 101–136, 1979. a, b
Wendland, H.: Scattered data approximation, Cambridge university press, vol. 17, https://doi.org/10.1017/CBO9780511617539, 2004. a
Wieczorech, M. A. and Meschede, M.: SHTools – Tools for working with spherical harmonics, Geochem. Geophy. Geosy., 19, 2574–2592, https://doi.org/10.1029/2018GC007529, 2018. a
Wimmers, A. J. and Velden, C. S.: Seamless Advective Blending of Total Precipitable Water Retrievals from PolarOrbiting Satellites, J. Appl. Meteorol. Clim., 50, 1024–1036, https://doi.org/10.1175/2010JAMC2589.1, 2011. a, b
Zender, C. S.: Analysis of selfdescribing gridded geoscience data with netCDF Operators (NCO), Environ. Modell. Softw., 23, 1338–1342, https://doi.org/10.1016/j.envsoft.2008.03.004, 2008. a
Zerroukat, M., Wood, N., and Staniforth, A.: A monotonic and positive–definite filter for a SemiLagrangian Inherently Conserving and Efficient (SLICE) scheme, Q. J. Roy. Meteor. Soc., 131, 2923–2936, 2005. a
Zerroukat, M., Wood, N., and Staniforth, A.: The parabolic spline method (PSM) for conservative transport problems, Int. J. Numer. Meth. Fl., 51, 1297–1318, 2006. a
Zienkiewicz, O. C. and Zhu, J. Z.: The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique, Int. J. Numer. Meth. Eng., 33, 1331–1364, 1992. a, b
 Abstract
 Copyright statement
 Introduction
 Remapping algorithms
 Metrics to evaluate remapping algorithms
 Metrics for Intercomparison of Remapping Algorithms (MIRA) workflow
 Results and discussion
 Future research directions
 Conclusions
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Copyright statement
 Introduction
 Remapping algorithms
 Metrics to evaluate remapping algorithms
 Metrics for Intercomparison of Remapping Algorithms (MIRA) workflow
 Results and discussion
 Future research directions
 Conclusions
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References