A system of conservative regridding for ice–atmosphere coupling in a General Circulation Model (GCM)

. The method of elevation classes, in which the ice surface model is run at multiple elevations within each grid cell, has proven to be a useful way for a low-resolution atmosphere inside a general circulation model (GCM) to produce high-resolution downscaled surface mass balance ﬁelds for use in one-way studies coupling atmospheres and ice ﬂow models. Past uses of elevation classes have failed to conserve mass and energy because the transformation used to regrid to the atmosphere was inconsistent with the transformation used to downscale to the ice model. This would cause problems for two-way coupling. A strategy that resolves this conservation issue has been designed and is presented here. The approach identiﬁes three grids between which data must be regridded and ﬁve transformations between those grids required by a typical coupled atmosphere–ice ﬂow model. This paper develops a theoretical framework for the problem and shows how each of these transformations may be achieved in a consistent, conservative manner. These transformations are implemented in Glint2, a library used to couple atmosphere models with ice models. Source code and documentation are available for download. Confounding real-world issues are discussed, including the use of projections for ice modeling, how to handle dynamically changing ice geometry, and modiﬁcations required for ﬁnite element ice models.


Introduction
Many questions still surround the issues of how ice sheets respond to climate forcing and how those changes will affect sea level, regional and global climate.Recent observations have shown accelerating mass loss from the Greenland and Antarctic ice sheets (Vaughan et al., 2013), adding urgency to these questions.Although general circulation models (GCMs) are able to project changes in ice sheet surface mass balance, projection of changes in ice sheet mass due to ice dynamics with GCMs remains a challenge (Church et al., 2013).A number of climate modeling groups are addressing these deficiencies by adding dynamic ice flow effects to existing GCMs: GISS ModelE (Schmidt et al., 2006), CESM (Hurrell et al., 2013), HadGEM2 (Collins et al., 2011), etc.This is done by coupling the GCM with an existing ice flow model, such as Glimmer-CISM (Rutt et al., 2009), BISI-CLES (Cornford et al., 2013), ISSM (Larour et al., 2012), PISM (Bueler and Brown, 2009), etc.A full understanding of the long-term evolution of an ice sheet within a coupled climate system requires coupling with the ocean as well as atmosphere.Surface runoff, ocean cavity circulation and salinity gradient effects are all important.In this paper, we focus only on coupling with the atmosphere.
One can distinguish between one-way and two-way coupling.In one-way coupling, the GCM is used to develop surface mass balance (SMB) and temperature fields, which are then used to drive the ice flow model off-line.This process misses effects caused by feedbacks from the ice sheet to the rest of the Earth system: for example, decreased ice sheet Published by Copernicus Publications on behalf of the European Geosciences Union.R. Fischer et al.: Conservative regridding for ice-atmosphere coupling albedo (Qu and Hall, 2006) or lowered atmosphere orography (Ridley et al., 2005).Past studies with one-way coupling have yielded useful insight into the future of present-day ice sheets; examples include Huybrechts (1994), Greve (2000), Stone et al. (2010), Bindschadler et al. (2013), Lipscomb et al. (2013), Nowicki et al. (2013a), Nowicki et al. (2013b) and Goelzer et al. (2013).However, ice sheet feedbacks are expected to be increasingly important for simulations of the long-term evolution of ice sheets and the climate associated with them.This kind of feedback probably plays a significant role in many events in the paleorecord (Dansgaard-Oeschger events, Heinrich events, the Younger Dryas).
Two-way coupling strategies address these issues by allowing the atmosphere to be influenced by changes in the ice sheet elevation, extent and albedo over time.We distinguish between loose and tight two-way coupling.Loose twoway coupling involves running a series of GCM simulations with different ice sheet configurations, each based on the result of the previous.Each GCM run is a separate simulation, without continuity of mass or energy between runs.Studies with loose two-way coupling have yielded insight into future equilibrium states for ice sheets and climate (Ridley et al., 2005(Ridley et al., , 2010)).However, ice sheets change configuration in these runs without accompanying mass and energy fluxes required to make those changes happen.This is equivalent to applying an unknown impulse forcing to the ice sheet each time it is changed.Although the coupled ice sheet might eventually reach the correct equilibrium state, the transients involved in achieving that equilibrium will be suspect.Unfortunately, results relevant to human society all require an understanding of the transients.For successful simulation of transients, we turn to tight two-way coupling.It involves running the ice flow model, step by step, along with the rest of the GCM -while conserving mass and energy along the way.Attention to conservation is required, since the GCM is simulating a more nearly closed system that could run for a long time.
When one couples dynamic ice flow models with GCM atmospheres, two models operating on different grids and timescales must communicate: ice flow models operate at low frequency on a high-resolution grid with local projection, while GCM atmospheres operate at high frequency on low-resolution global grids.A number of issues arise due to this mismatch, including how one creates high-resolution surface mass balance fields from low-resolution GCM input.Elevation classes address the latter issue.They were first introduced for precipitation downscaling in a GCM by Leung and Ghan (1998) and later applied to one-way coupling from GCM atmosphere to ice flow models by Lipscomb et al. (2013).The key insight is that mass and energy fluxes between the atmosphere and an ice sheet vary approximately by elevation within a local region.
When using elevation classes, a third grid is introduced, the elevation grid.This allows the GCM to compute surface fields at a variety of elevations within each atmosphere grid cell, not just the elevation seen by the atmosphere.A highresolution surface mass balance is produced on the ice grid by first computing SMB on the elevation grid, and then using a vertical interpolation scheme to produce SMB on the ice grid.This method of interpolation produces surprisingly good results: although it cannot capture certain localized effects (e.g., wind direction), it has been shown to allow GCMs to produce surface mass balance fields approaching the quality of those produced by regional climate models (Vizcaino et al., 2013).
Inside the GCM, SMBs computed on the elevation grid must be regridded to the atmosphere grid as well as the ice grid.In order to maintain conservation in a tight two-way coupled system, it is essential that the set of regridding operation chosen is self-consistent: that is, if a flux field on the elevation grid is regridded simultaneously to the atmosphere and ice grid, then the total amount of flux represented by the resulting two fields should be the same.For conservation purposes, the specifics of these two transformations are not important, as long as they are consistent with each other.Past efforts at one-way coupling have defined these two transformations in ways that each make intuitive sense, but are not consistent with each other.This is not a problem for one-way coupling, but it would cause conservation problems for twoway coupling.
This paper develops the concept of the elevation grid on which the ice surface model runs, and then derives a set of conservative regridding transformations between the atmosphere, elevation and ice grids.The coupled processes under consideration along with the transformations required for tight two-way coupling are introduced in Sect. 2. Section 3 focuses on the elevation grid, while the grid fundamentals necessary for the two-way coupling are presented in Sect. 4. Section 5 deals with the use of projection and associated issues encountered when bridging between the spherical geometry of GCMs and Cartesian ice sheet models.We show how to choose and implement the transformations in Sects.6 through 9 and work through realistic examples of these transformations in Sects.10 through 12.We touch on a number of extra "wrinkles" in the real-world problem: procedures to use when the elevation grid is based on a horizontal grid other than the atmosphere grid in Sect.13, and regridding procedures required when ice elevations, ice extent or elevation grid change in the simulation in Sect.14.Finally, we present in Sect.15 a library, Glint2, that can be used to tightly couple GCMs and ice sheet models.

Coupled processes
The coupled atmosphere-ice system involves three models interacting with each other: an atmosphere model, an ice flow model and an ice surface model situated between them (Fig. 1).The atmosphere model is coupled with the ice surface model, which tracks the top few meters of ice.Processes modeled here include a full surface mass-energy balance computation, snow-firn compaction, albedo effects, meltwater percolation, runoff, refreezing, etc.The bottom of the ice surface model is coupled to the ice flow model.

Ice surface model
In order to couple an ice sheet to a GCM atmosphere, SMB on the ice sheet must be calculated from GCM outputs.Some one-way coupled studies have used temperature index schemes (Huybrechts, 1994;Greve, 2000;Stone et al., 2010;Bindschadler et al., 2013;Nowicki et al., 2013a, b;Goelzer et al., 2013): the mean surface temperature over each coupling time step (typically one month or year) is used to compute SMB, and both are passed to the ice flow model.It is hard to see how energy can be conserved with such a scheme.The problem is that the atmosphere must run for many time steps before atmosphere-ice sheet energy fluxes are computed on a coupling time step.The computed flux on the ice sheet will not be the same as that sum of fluxes seen by the atmosphere over the previous coupling time step -and it is "too late" to go back and change the atmosphere to match.
For this reason, a full energy balance scheme, computed each atmosphere time step (typically one hour) and integrated over the coupling time step, is considered essential in a GCM setting.Energy flux between atmosphere and ice sheet follows diurnal and seasonal cycles, making positive and negative contributions to the integrated flux.It is important that energy flux is computed at a small enough time step to capture these effects accurately; typically, the atmosphere time step is sufficient.
While ice sheets move slowly and can be modeled with long time steps, the modeling of the surface of the ice sheet requires short time steps.This is accomplished by introducing the ice surface model as a third model, sitting in between the ice flow model and atmosphere model.The top of the ice surface model couples with the atmosphere model every atmosphere time step, whereas the bottom couples with the ice sheet model every coupling time step.The ice surface model needs to be thick enough so there is little variation in temperature at the bottom: 15 m is sufficient.This ensures that the heat flux with the ice flow model will be small and low frequency.
The ice surface model can serve an additional purpose of modeling the great variety of surface processes that may be relevant to the long-term evolution of ice sheets: snow-firn compaction, surface runoff/drainage networks, water percolation, refreeze, albedo effects and wind-blown snow, to name a few.
The use of an ice surface model solves some important problems, but it also introduces nonphysical elements into the model.Ice contained in the ice surface model does not advect along with the ice flow model, nor does it contribute to the stress field of the ice below it in the ice flow model.In both cases, we expect the relative error to be small: the top 15 m of snow/firn contains less than 1 % of the mass of a 1000 m thick ice sheet.There may be ways to fix these problems -however, there is no need to make the ice surface model thicker than it needs to be.Numerous studies have shown that ice sheets below about 15 m are fully insulated from surface weather and seasonal cycles (Zagorodnov et al., 2006).

Three models, three grids
Each of the three models in the coupled system runs on its own grid.The atmosphere is run on the atmosphere grid (A) and the ice flow model is run on the ice grid (I ).The ice surface model is run on the elevation grid (E), which is based on the atmosphere grid (see Sect. 3).All three grids are twodimensional, in the sense that they are used to construct twodimensional functions f (x, y) over the domain.Regridding operations are needed to pass mass and energy fluxes between the models.
It is important to keep in mind the relative size of the three grids.We set up a test using the GISS 2 • × 2 1 2 • atmosphere grid (Schmidt et al., 2006), overlapping with the 5 km grid from SeaRISE (Sea-level Response to Ice Sheet Evolution; Bindschadler et al., 2013).We used 40 elevation points, spaced every 100 m from 0 m to 4000 m.In this case, A had 146 grid cells, I had 66 906 and E had 1829.These numbers only account for grid cells involved with the Greenland ice sheet.In general, the ice grid will be finest, the atmosphere grid coarsest and the elevation grid in the middle.Time frequency mismatches are another issue.The atmosphere runs at high frequency, each time step being typically 1 h.On the other hand, the atmosphere and ice flow models are coupled at much lower frequency, typically one month or any other time period.We call these two time steps the atmosphere time step and ice coupling time step.

Fully coupled system
Figure 2 shows the data flow of the fully coupled system.Steps of the data flow are organized based on their frequency: the top circle of steps runs at the same frequency as the GCM atmosphere, while the bottom circle runs at the ice coupling frequency -typically one month or more.We describe the steps involved in coupling the GISS ModelE (Schmidt et al., 2006) with PISM (Bueler and Brown, 2009); however, these steps are general for any GCM or ice flow model.We now trace through the data flow on a typical coupled run.

Atmosphere time step
When the atmosphere runs, it produces a set of fields on the atmosphere grid that affect the processes in the ice surface model: for example, downwelling long-wave radiation, downwelling shortwave radiation, precipitation, etc.These fields must be regridded from A to E. We denote the set of fields to be regridded with a capital letter, R A ; the result of the regridding is R E (see Appendix A and Fig. 3 for notational conventions).
The ice surface model is run at the same frequency as the atmosphere.Among its outputs are mass and energy fluxes with the atmosphere: evaporation, sublimation, upwelling long-wave radiation, latent heat release, etc.These are represented by −F E and are regridded to −F A before being passed back to the atmosphere on the next time step.
The ice surface model also produces fluxes in the ice flow model; these are described immediately below.

Coupling time step
On each atmosphere time step, relevant flux outputs from the ice surface model are accumulated as FE for future coupling with the ice flow model.They are named FE because these fluxes are in general equal and opposite to fluxes sent to the atmosphere.Every coupling time step -about once a month -the accumulated FE is regridded to the ice grid ( FI ) and passed to the ice flow model.
The ice flow model produces changes in ice surface topography and extent, as well as a small energy flux between the ice flow and ice surface models (together, we call these D I ).Changes in ice topography and extent are regridded to the atmosphere grid (D A ) and used to adjust the atmosphere's orography.The energy flux is regridded to the elevation grid (D E ) and applied to the ice surface model.

Regridding requirements
We see that the fully coupled system requires the use of five regridding operations at various points in its run (Fig. 2).These regridding operations must be conservative, in the sense that none of them can change the integral of the field on the domain as a whole.In fact, we would like to impose a stronger conservation condition so that values are conserved within each atmosphere grid cell.
We develop these five regridding operations in the sections below.We begin by discussing the method of elevation classes in a general manner (Sect.3) and then move on to developing basics of conservative regridding (Sect.4).Because conservation is defined in terms of integration over areas, part of our discussion involves a definition of how to integrate functions on the elevation grid E. Finally, we show how the required set of operations can be constructed in a fully self-consistent manner.
Throughout, we assume that the atmosphere grid A has L0 parameterization: functions on A are represented by their mean value within each grid cell.This becomes our basis for conservation: all regridding operations are conservative over every atmosphere grid cell.We do not require any specific parameterization for the ice grid, but we consider cases for L0 (piecewise constant) and L1 (piecewise linear) parameterizations, the latter which are commonly used in finite element ice models.

Elevation points
The method of elevation classes, when applied to the atmosphere-ice coupling problem, involves running the ice surface model at one or more elevations in each atmosphere grid cell (Lipscomb et al., 2013).Temperature, pressure and precipitation are extrapolated to a set of elevations within the grid cell.They are then used in an ice surface model to compute a full surface mass and energy balance at each elevation.The result is a set of "what-if" scenarios, giving an estimate of what the fluxes would have been, had the ice surface of the grid cell been at the given elevation -rather than the elevation seen by the atmosphere model.
The modeler must choose which elevations to use for each atmosphere grid cell.The simplest approach is to use a fixed set of elevations across all grid cells -for example, every 100 m from 0 m to 4000 m.
However elevations are chosen, the result is a new "grid" -the elevation grid -on which the ice surface model is run and surface fluxes are generated.The elevation grid is derived from the atmosphere grid, in the sense that each elevation grid "cell" (or elevation point) is associated with one parent atmosphere grid cell.If elevation point E j is associated with atmosphere grid cell A i , we write E j ∈ A i .
Once the GCM has computed a conserved quantity on the elevation grid, those values can be used to develop a relation between elevation and SMB within each atmosphere grid cell.Suppose we have computed a flux field f E : SMB, for example.For an atmosphere grid cell A i , consider the components of the flux field f E j that are related to the enclosing atmosphere grid cell A i .Or more formally, consider f E j for all j such that E j ∈ A i .We can think of these component  values as samples of a 1-D function relating elevation to flux (SMB), within the localized region of A i .By interpolating between those points using standard methods, one can construct a continuous function relating flux to elevation within A i .As long as enough points are used and the function being interpolated is smooth enough, this procedure will yield an arbitrarily accurate representation of the "true" function.
In fact, the functions we typically wish to interpolate -surface mass and energy balance averaged over about a monthare quite smooth as functions of elevation (Fig. 4).In the face of spatially invariant precipitation, one would expect SMB to be constant to first order above the equilibrium line altitude (ELA), and to decrease linearly with elevation below the ELA.Near the ELA, one would expect a smooth transition because the ELA goes up and down over a month of diurnal cycles.This is in agreement with experimental work, which has shown SMB below the ELA to vary linearly with temperature (Braithwaite, 1981;Box et al., 2013).
How many elevation points are required to properly resolve the elevation-flux relationship?This depends on the interpolation scheme: higher order schemes will require fewer points than piecewise linear interpolation.But the general shape of the function -two straight lines connected by a curve near the ELA -implies that not many points are needed, except for near the ELA.
If one knows where the ELA is on every atmosphere grid cell, then this is a useful guide in selecting elevation points.But if one is studying ice sheets in a changing climate, then the ELA will be expected to move over time.It is possible to adaptively move elevation points as the ELA moves.But it is simpler just to use many points everywhere.In our tests, we have used 40 points at 100 m spacing, which is probably more than sufficient for coupled ice sheet simulations.

R. Fischer et al.: Conservative regridding for ice-atmosphere coupling
We have not done a careful study of the optimal number of elevation points.
The method of elevation classes -or elevation pointsprovides a way to construct an interpolated relation between elevation and surface mass-energy flux within each atmosphere grid cell.Additional choices need to be made in order to produce fully downscaled flux fields on the ice grid (Sect.6).We are now almost ready to define the transformations posited in Sect.3.But first, we pause for some groundwork on grid fundamentals in Sect. 4 and projection issues in Sect. 5.

Grid fundamentals
Numerical models represent continuous fields with linear combinations of a finite set G of basis functions -which we call a parameterization (Appendix A describes our mathematical notation).For example, suppose that G uses the basis functions g(x, y) = [g 1 (x, y), . . ., g n (x, y)], where g(x, y) is the vector of all the basis functions.Suppose we have an n-dimensional vector f G with components f G i .That vector represents the function f G (x, y) formed by a weighted sum of the basis functions: A wide variety of basis function sets are used for different problems.Most commonly in climate modeling, f G i represents the mean value of f G (x, y) within some well-defined bounded region, which we call a grid cell.From a conservation point of view, a function f G (x, y) with L0 parameterization may be taken to be constant within each grid cell, with discrete "jumps" from one grid cell to the next.In this case, the basis function g i (x, y) for grid cell i is equal to 1 inside the cell and 0 outside.L0 parameterizations are widely used in climate and finite difference ice flow models.Finite element ice flow models might use L1 or even higher order parameterizations (Zienkiewicz et al., 2013), and they use the term "mesh" to describe the geometry of their basis functions.In this paper, we use the term "grid" to refer to both the vector space in which f G lives and its associated parameterization.

Integration on grids
Because we are working with conserved quantities, it is essential that we can integrate functions over any well-defined region B. Because a parameterization defines f G (x, y) on every point, this integral is well-defined.Substituting from Eq. (1) and using linearity, we get If we wish to integrate over B repeatedly, we can precompute B g(x, y)dA, and then take dot products with f G as needed.How one integrates the basis functions depends on the nature of those basis functions.We give specifics for L0 and L1 grids in Appendices B and C, respectively.

Comparison across grids
In regridding applications, we need to compare fields across different grids.They cannot be tested for simple point-bypoint equality, due to differences in grid structure.Instead, we compare by integrating two fields over the same region or set of regions.
If we have two fields f G and f H on two different grids, we say that Suppose we wish to compare f G and f H over an entire domain?If we have a set of regions A = {A 1 , . . ., A n } tiling that domain, then we can say the two fields are equivalent on A if they are equivalent on all of A 1 , . . ., A n .Note that A could be an L0 grid, or simply a set of regions on the domain.If we have a regridding transformation F()

Area-weighted remapping
Suppose we have two grids G and H , with H being L0 parameterized -for example, G might be an ice grid and H an atmosphere grid.Suppose we have a field f G , which we wish to regrid in a conservative manner to result in the field f H .One common way to do this is to compute each component f H i based on the value of f G (x, y) integrated over the area covered by the ith grid cell in the grid H . Or, more formally, By construction, this transformation is conservative on H .Note that f G and f H will generally not be equivalent over other regions other than grid cells in H -for example, over grid cells in G (if G happens to be L0 parameterized).
Area-weighted remapping is closely related to our definition of equivalence above.If we have two fields on two different grids, equivalence on H can be tested by regridding both fields to H and comparing the resulting vectors.The transformation from G to H is linear and thus representable by a matrix.It is variously called area-weighted remapping or conservative regridding (Ramshaw, 1985).The matrix is computed by integrating the source basis functions g 1 (x, y), . . ., g n (x, y) over every grid cell in the destination grid.

Atmosphere
Grid Ice Grid Exchange Grid Fig. 5.An exchange grid, obtained by overlapping a sample atmosphere grid and ice grid (not to scale).Each resulting grid cell, irregularly shaped, overlaps with exactly one atmosphere grid cell and exactly one ice grid cell.

Interpolation grid
We are now almost ready to discuss our implementations of the transformations E → I and E → A (see Appendix A and Fig. 3 for mathematical notation).But first, we must introduce the interpolation grid G, which is used to rigorously define basis functions on the elevation grid E. The user has a choice of how G is to be chosen.
In general, G is chosen simply as the ice grid I .If I is L0-parameterized, we might instead choose G to be the exchange grid between A and I (Fig. 5): the L0 grid whose grid cell outlines are formed by the intersection of grid cells in A and I (Balaji et al., 2006).The exchange grid is a useful choice for G because every exchange grid cell overlaps at most one atmosphere grid cell.This choice has its pros and cons, which we explore in Sect.9. Either way, the interpolation grid G is similar to the ice grid I and can be thought of as an ice grid proxy in most cases.

Projection issues
Whether the source grid is L0 or L1, area-weighted remapping algorithms need to find the intersection of grid cell outlines from two grids.Technically, this is only possible if the two grids exist on the same surface.In our case, atmosphere models exist on the surface of a sphere, whereas ice flow models work on a Cartesian plane.Unless an ice flow model formulated in spherical coordinates is used, the exchange grid between an atmosphere grid and an ice grid cannot be directly computed.
This problem is solved using a map projection (Snyder, 1987).We let A be the original atmosphere grid and then let A be the projected atmosphere grid -projected to the Cartesian plane using the chosen projection.Regridding computations described in this paper are made to/from A, the projected grid.
In general, the area of a grid cell can change when it is projected.Changes in area due to the use of a projection that does not preserve areas are called projection error (Lauritzen and Nair, 2008).Even if an equal area projection is used, in practice grid cells still experience small changes in area.This is because projected grid cells have complex shapes, but the algorithms in Appendix D are only able to approximate these shapes with polygons.Changes in area due to polygonal approximation of curves, or numerical artifacts in the methods used, are called geometric error (Lauritzen and Nair, 2008).
Whether a change in grid cell area arises from projection or geometric error, data must be rescaled to maintain conservation when transforming a field f A to f A : (5) Note that this transformation is diagonal and invertible.Because this is simply a rescaling, the regridding schemes presented in this paper are not operationally affected by projection issues.
Although this rescaling scheme may be used to address any change in grid cell area, it is preferable to eliminate or minimize these changes to begin with.To this end, we consider projection and geometric error separately.

Projection error
The stereographic projection used in SeaRISE (Bindschadler et al., 2013) can change the area of grid cells by up to 10 %.If a projected grid cell is 10 % smaller then the original, then that means that 1 m of accumulation in the GCM will turn into 1.1 m in the ice flow model.This could cause significant discrepancies in dynamic ice flow.Because the projection in SeaRISE is constructed to minimize errors in a band along the 71 • N parallel, we would not expect the magnitude of projection error to change significantly if an oblique stereographic projection were used instead.
One solution to the problem of projection error would be to use an equal area projection -the Lambert equal area projection, for example.However, equal area projections do not preserve angles, which distorts the ice dynamics at a local scale.We conclude that nonphysical distortions happen whether an area-preserving or angle-preserving projection is used.It is not yet clear which choice gives better results in the end.
One innovative approach to this problem is to use an anglepreserving projection, but to allow ice flow model grid cells to vary in size, based on the local distortion of space introduced by the projection.Parameters dx and dy are kept for each grid cell.These grid parameters are included in the ice flow model equations, thereby accounting for shape and area distortion caused by the projection in a physically meaningful way (Pollard and DeConto, 2012).

Geometric error
Geometric errors are typically three orders of magnitude smaller than projection errors.In theory, they may be Fig. 6.An example of geometric error, in which grid cells change size due to polygonal approximation.This map, made using a Lambert equal area projection centered on the North Pole, shows a set of latitude-longitude grid cells -a kind commonly used in atmosphere models.Solid blue lines show polygonal approximations, while dotted red shows the actual grid cells on the sphere.Note that all grid cells (in this local map) shrink when approximated.
arbitrarily reduced further by increasing the number of sides used in the approximating polygons -the parameter n specifies the number of line segments used to approximate each side of the original grid cell in its projected state.But this method has its limit in practical terms.
Consider a typical latitude-longitude grid on the sphere with a Lambert equal area projection centered at the North Pole (Fig. 6).In this case, the area of a polygonal approximation will always be smaller than the area of the actual grid cell: we end up inscribing polygons inside of circles.The geometric error will depend on the number of sides of the inscribed polygon, which in this case depends on the number of grid cells in the circle, and on n.This is the method used by Archimedes to approximate the value of π in ∼ 250 BCE (Heath and Archimedes, 1897, p. 91).Unfortunately, it converges only quadratically, as O(1/n 2 ).Meanwhile, memory use to store all those polygons goes up by O(n).Memory and time requirements will therefore be exponential in terms of the number of digits of accuracy required: each additional digit of accuracy will require an increase in the number of sides by a factor of √ 10.It is therefore not practical to make geometric error arbitrarily small by increasing n.In our experience, a value of n = 2 yields geometric error of approximately 10 −4 .We believe n = 2 offers an acceptable trade-off between efficiency and accuracy.

Transformation: elevation to ice/interpolation grid
We can now derive the first of the five transformations posited in Fig. 2. We begin with the transformation E → G from the elevation grid to interpolation grid, because this is the transformation that ice modelers are most concerned with when seeking high-quality downscaled flux fields.In various contexts, this transformation might be referred to as an interpolation, downscaling or regridding operation.
In Sect.3, we showed how values f E at elevation points may be used to construct a flux-elevation relationship f i (z) in each atmosphere grid cell A i .These relationships may then be used to construct f G on the interpolation grid.The resulting transformation for E → G potentially involves horizontal and vertical interpolation, of which the modeler has considerable choice.
We present here three methods of interpolation from E to G: z interpolation (Sect.6.1), bilinear interpolation (Sect.6.2) and elevation class interpolation (Sect.6.3).All three methods assume that G is L0 parameterized and can be used for the basis of a conservative coupling.And since they are all linear, they can all be represented by a (sparse) matrix, which we will call M. The methods are extended to the L1 case in Sect.6.4

Z interpolation
Z interpolation constructs a flux field f G on the interpolation grid through direct application of the flux-elevation relationship f i (z) derived in Sect. 3 for atmosphere grid cell A i .Supposing G is L0 parameterized, consider an interpolation grid cell G j , wholly contained in one atmosphere grid cell A i , with mean elevation z j .We would set the value of that grid cell to f G j = f i (z j ).If G j intersects more than one atmosphere grid cell, the same procedure is followed for each atmosphere grid cell A i it intersects -and the results are summed together, area-weighted by |G j ∩ A i |.
We call this z interpolation.As long as standard interpolation methods are used to construct f i (z), z interpolation defines a linear function from f E to f G and can be represented by a matrix.An example of the result is shown in Fig. 7a.In general, z interpolation produces SMB fields that vary smoothly within each atmosphere grid cell but that contain discontinuities between grid cells.

Bilinear interpolation
There is some concern that discontinuities created by z interpolation could cause problems as an input to an ice flow model.In that case, a bilinear interpolation step may be used to create a smooth field, as in Lipscomb et al. (2013).This scheme sets the value of interpolation grid cell G j equal to a linear combination of f i (z), f k (z), f l (z) and f m (z), where atmosphere grid cells A i , A k , A l and A m are the four grid cells with centers closest to the center of G j (Fig. 8).Results are shown in Fig. 7b.Bilinear interpolation has an advantage over z interpolation in that it produces smooth fields.However, bilinear interpolation presents a number of problems: -The fields it produces will have a significantly different total mass than the fields produced by z interpolation.
Our experience with monthly SMB fields over Greenland indicates that storms tend to leave large amounts of snow in a few localized areas.Bilinear interpolation tends to reduce the total amount of snowfall in these cases, causing potentially significant differences in the GCM model run.
Variations in total mass caused by the choice of interpolation procedure will not cause conservation problems: we make this clear in Sect.7.However, an interpolation scheme that produces a significantly different mass from the apparent "intent" of the GCM could cause biases.In particular, a scheme that makes snowstorms smaller than the GCM "intended" would produce a negative bias in the equilibrium extent of the ice sheet.
-It introduces significant numerical diffusion into the system.
-The A → E transformation derived with bilinear interpolation can introduce nonphysical artifacts (Sect.12).-It is also not always clear how to extend bilinear interpolation to the case of non-rectangular atmosphere grids.This problem can also affect non-regular points in mostly regular grids (e.g., at the poles of a latitudelongitude grid).

Elevation class interpolation
One final choice for regridding is to define elevation classes as they were originally formulated (Leung and Ghan, 1998).
In this approach, each atmosphere grid cell is grouped into sub-regions based on elevation bands.For example, an atmosphere grid cell near the coast spanning elevations of 0 m-600 m might be grouped into sub-regions of 0 m-200 m, 200 m-400 m and 400 m-600 m.The grid E can then be thought of as having irregularly shaped grid cells formed by the intersection of atmosphere grid cell boundaries and elevation contours (Fig. 9).Area-weighted remapping is used to interpolate from the elevation grid to the atmosphere grid.Elevation class interpolation is equivalent to a specialized form of z interpolation, in which the 1-D function f i (z) (for grid cell A i ) is interpolated as a discontinuous piecewise constant function (Fig. 10), rather than a piecewise polynomial (Fig. 4).Unless one expects significant discontinuities in topography or the elevation-flux relationship, such a zero-order interpolation scheme would be expected to give a less accurate version of f G than a higher order scheme such as z interpolation.Since we expect the elevation-flux relationship to be smooth with elevation, we recommend the use of z interpolation over elevation class interpolation.

Interpolating to L1 grids
We have outlined methods for the interpolation E → G, as long as G uses an L0 parameterization.The procedures need to be modified for interpolation or ice grids using L1 parameterization.In a finite element mesh ("grid"), field values are determined at mesh vertices and linearly interpolated within each triangular element (Zienkiewicz et al., 2013).Any of the interpolation methods above may be used to determine the value of f G (x, y) at each vertex.Once vertex values have been interpolated on the vertices of a finite element mesh, the value of f G (x, y) at all other points is fully determined.

Transformation: elevation to atmosphere grid
As shown in Fig. 2, a regrid step from the elevation to the atmosphere grid is required on every GCM time step.Once the modeler has chosen the transformation E → G, we show here how to derive a transformation for E → A that is consistent with it.
To derive this transformation, consider one GCM time step, during which a flux field f E between the atmosphere and ice surface is computed on the elevation grid.That field will be regridded both to the interpolation and atmosphere grids.Conservation requires that the resulting fields are equivalent on A, i.e, If the transformation E → G derived in Sect.6 is represented by the matrix M, then we can write f G = Mf E .Substituting into Eq.( 6), we get the requirement Equivalence on A can be tested by conservatively remapping both sides to A and then comparing (in this case, the left-hand side is already on A, so no remapping is required there).Denoting R as the matrix that conservatively regrids from G to A, we have the requirement If we treat this as a definition for the transformation E → A, we have derived a transformation for E → A that is consistent with the transformation we already chose for E → G.
In other words (Fig. 11): we will regrid from E to A by first regridding from E to G (represented by the matrix M) and then use area-weighted remapping from G to A (represented by the matrix R).The GCM, which needs to use the transformation E → A, does not need to know anything about the ice grid used to derive that transformation.It

Ice Grid I Elevation
Grid E

Area-Weighted Remapping
Atmosphere Grid A Fig. 11.The transformation from elevation grid to atmosphere grid is derived by first interpolating to the ice grid (represented by matrix M) and then using a conservative area-weighted remapping step to the atmosphere grid (matrix R).The two steps may be combined by computing the matrix product RM.This construction ensures that a quantity computed on the elevation grid will have the same total mass when regridded to the ice grid or atmosphere grid.The two transformations from the elevation grid to the ice and atmosphere grids are said to be consistent.
just needs to know the final matrix RM, which can be precomputed via matrix multiplication.
We have derived a transformation for E → A that by definition is consistent with a previously chosen transformation for E → G. Our approach differs from previous efforts, which would start out with E → A and try to find a transformation for E → G that is consistent with it.

A basis for the elevation grid
Every vector space has basis functions, including the elevation grid E. Our choice for the transformation E → A constrains the basis functions we use on E. To find these basis functions, we expand on the central idea in the previous section, i.e., that we can determine properties of a field on the elevation grid by regridding to the ice grid.In this case, we define f E (x, y) to be equal to f G (x, y), where f G = Mf E .In other words, we can evaluate f E (x, y) at a point by interpolating f E to G and then evaluating f G (x, y).This definition is consistent with our method for f A in the previous section.
We now use this principle to obtain a formula for the basis functions of E. Substituting into Eq.( 1), we get Rewriting in indicial notation, using the commutative property of multiplication, and swapping index letters, we get We can use this, along with Eq. ( 1), to determine the basis functions for the elevation grid: or, switching back to vector notation, e(x, y) = M T g(x, y).
Equivalently, the basis function e i (x, y) is the function obtained if we set the f E i = 1 and all other components of f E = 0, regrid to the interpolation grid, and then examine the resulting function f G (x, y).
In general, these basis functions will be not orthogonal.Their exact form depends on choices made in choosing M, i.e., vertical and horizontal interpolation choices, as well as grid geometry issues.We have plotted some example basis functions in Fig. 12.
Note that if elevation class interpolation is used (Sect.6.3) and G is the exchange grid, then our scheme reduces to a traditional elevation class scheme, and basis functions will represent constant-value sub-grid tiles, which are orthogonal.

Forward transformations: summary
We have now defined three of the five regridding transformations required by Fig. 2: E → G, G → A and E → A. We call these the forward transformations because they are linear and can be represented by matrices.See Fig. 13 for a diagram of how these transformations may be used to regrid fields.
We defined these three transformations in a consistent manner.We began by allowing the user choice in constructing the transformation from the elevation to interpolation grids, represented by matrix M: this makes sense because users desire choice for E → G.We then noted that standard area-weighted remapping, represented by matrix R, may be used to regrid from the interpolation grid to the atmosphere grid.From these two choices, transitivity requires that the transformation from the elevation to atmosphere grid is represented by the matrix RM.Conservation is maintained because E → G and E → A are consistent with each other, producing f G and f A that are equivalent on every grid cell in A. Our choices for these transformations lead directly to a well-defined set of basis functions for the grid E, examples of which we plotted.Fields on this grid will be represented in terms of these basis functions -for example, surfaceatmosphere fluxes and ice surface model state.

Reversing the transformations
We have defined three transformations so far, but Fig. 2 indicates that five are necessary for full functioning of the coupled system.We still need to derive appropriate procedures for the "reverse transformations" A → E and G → E, indicated in Fig. 13 as dotted lines.Because of the differences in dimensionality between the three grids (Sect.2.3), M and RM are not invertible in a simple linear algebraic sense: A → E is underdetermined and G → E is overdetermined.We must still find ways to compute these transformations that retain conservation over atmosphere grid cells.

Transformation: atmosphere to elevation grid
Suppose we have a flux field f A on the atmosphere grid A.
We wish to regrid it to an "equivalent" flux field f E on the elevation grid E. In this case, f might represent precipitation or downwelling radiation from the atmosphere.
The problem is underdetermined.Any solution for which f A ≡ A f E will be conservative.This is the same as requiring that Because of the many-to-one relationship between elevation and atmosphere grid cells, our intuition tells us that f E may be constructed simply by repeating values of f A within each atmosphere grid cell.This is physically selfconsistent.More precisely, if E j ∈ A i , then we would like to set f E j = f A i .We define a simple linear transformation that does this: f E = f A .Surprisingly, this definition for A → E is only conservative in some cases.If RM computes each f A i as a weighted sum only of values f E j where E j ∈ A i , then we say that RM is a local transformation -it uses only data from "within" an atmosphere grid cell to compute a value on that grid cell.In that case, it is easy to prove conservation by showing that RM f A = f A (see Eq. 13).The weights involved in RM do not matter, as long as they sum to 1 for each component of f A .
Therefore, if RM is local, we can simply use for our transformation A → E.Even better, it is easy to show that there will be no numerical dispersion in the round-trip transformation E → A → E because RM is the identity matrix.The lack of dispersion on this round-trip is important because it is computed every atmosphere time step.

Nonlocal RM
If RM is not local, then in general f A and f A will not be equivalent on A or even the entire ice sheet.We are forced to trade off between the most physically self-consistent value for f E = f A vs. something that conserves.We can use , along with quadratic optimization, to guide us to such a compromise.
We seek a vector f E such that We also wish that vector to be as close as possible to our intuition above.That is, we wish to minimize the quantity where our L 2 norm is weighted by the weight of each basis function.That is, we use a weight vector w where This is a sparse quadratic optimization problem with equality constraints.A number of numerical packages can solve it; we used GALAHAD (Gould et al., 2003).In our tests, solution typically requires a fraction of a second on a single core.This is so fast that we have found no need to consider suboptimal solutions to this problem.This procedure introduces some numerical dispersion into the fully coupled system from nonlocal regridding operations; by numerical dispersion, we mean movement of mass between adjacent atmosphere grid cells, thereby violating our desired property of conservation with each atmosphere grid cell.Considering the round-trip transformation E → A → E, it is clear that this transformation is not local, both because RM is not local, and because the quadratic program set up for A → E will be nonlocal.

Practical issues for nonlocal RM
GCMs are well equipped to deal with sub-grid tiles, each one occupying a fraction of the overall grid cell.The GCM will typically implement E → A by computing a weighted sum over sub-grid tiles in each grid cell, with weights based on each tile's fractional area.This capability is used to implement traditional elevation class schemes.
If the RM matrix is local, then by definition the value on each atmosphere grid cell is a weighted sum of the elevation points in that grid cell.This is compatible with existing GCM practice that assumes sub-grid tiles.To use the methods in this paper, the computation of E → A inside the GCM does not need to be replaced, it only needs to be fed a new set of weights.Even though elevation points do not have well-defined areas, the weight for elevation point j can be thought of as the "fractional area" that an elevation point contributes to its containing atmosphere grid cell.The GCM can be coded as if traditional elevation classes were being used, even if the user has chosen a more numerically accurate form of vertical interpolation in the construction of the M matrix.
Things get more complicated for the GCM if the RM matrix is not local.Instead of using a simple set of weights, the GCM will have to be multiplied by a general sparse matrix RM when computing E → A. An MPI (Message Passing Interface) gather will be required to compute A → E. Finally, some GCMs use implicit schemes for ice surfaceatmosphere coupling (Best et al., 2004), requiring a matrix inversion on every time step.If RM is nonlocal, then a large, sparse, global matrix inversion will be required, rather than a number of small, local inversions.Because of the significant practical complications that arise with the use of a nonlocal RM matrix, most users find it simplest to use a local RM matrix if at all possible.

When is RM not local?
The use of a nonlocal RM transformation can be an added burden to the GCM developer.It is therefore important to delineate cases in which RM is not local.
For L1 ice grids, as used with a finite element ice flow model, RM will be nonlocal: ice elements that straddle two atmosphere grid cells will by necessity involve elevation point values from two atmosphere grid cells.This will make R slightly nonlocal.
Even for L0 grids, M could be nonlocal, depending on the choice of interpolation schemes for E → G. Bilinear interpolation is inherently nonlocal, whereas the other interpolation strategies mentioned above are local.
For L0 ice grids with a local interpolation scheme, RM can be made local, as long as G was chosen to be the exchange grid (see Sect. 4.4).In this case, M will be local -because each exchange grid cell attains its value from elevation points in just one atmosphere grid cell.R will be local as well, for the same reason.We wish to regrid it to an "equivalent" flux field f E on the elevation grid E. As before, we will set up a quadratic optimization problem.
It would be nice if we could find f E such that Mf E = f G .This will not usually be possible, since G has far more degrees of freedom than E. Instead, we will minimize the quantity while we maintain conservation on each atmosphere grid cell (where R defines area-weighted remapping from G to A; see Fig. 13): This quadratic optimization problem may be solved with the same methods as in Sect.8.2.We weight components of our L 2 norm by the integral of each basis function in G.

Exchange grid or ice grid?
So far, we have defined our transformations in terms of the interpolation grid G -which could be either the ice grid or the exchange grid.But real ice flow models operate on the ice grid, and transformations must ultimately transform to/from that grid.If we have chosen G as the exchange grid, we need to extend our transformations above to regrid to/from the ice grid.We do this by using area-weighted remapping as necessary between G and I .
More formally, we construct transformations for X : G → I and X : I → G using area-weighted remapping.We can then represent E → I as the matrix product XM.Similarly, we can construct I → E by first computing f G = X f I and then using the reverse transformation G → E to obtain f E (see Fig. 13).
Note that the transformation X : G → I is not conservative on A. For this reason, the transformation E → I represented by XM is not conservative on A either, although it is conservative overall: "mass" lost from one atmosphere grid cell will be gained by neighboring cells.This will cause numerical dispersion when these transformations are used in the fully coupled system (Fig. 2).
Assuming the ice grid is L0 parameterized, we will now address the practical issue faced by the user, namely whether to choose the exchange grid or ice grid as G.If one uses the exchange grid for G in conjunction with a local vertical interpolation scheme, then the RM matrix will be local.This has a number of advantages.It eases implementation (Sect.8.3).There will be no numerical dispersion for the transformation E → A, used every atmosphere time step.The reverse transformation A → E, also used every atmosphere time step, will be trivial.However, there will be numerical dispersion for the transformation E → I , which is used once every ice coupling time step.
If one directly uses the ice grid for G, then the RM matrix will not be local, producing numerical dispersion for the transformation E → A and complicating the reverse transformation A → E. However, there will be no numerical dispersion for the transformation E → I .
On balance, we recommend keeping the RM matrix local if possible.Not only does this simplify implementation, it moves unavoidable numerical dispersion away from the atmosphere time step to the less frequent ice coupling time step.

Exchange grids for finite elements
The discussion in Sect.9 assumes an L0 ice grid.The results can be extended to L1 or higher order parameterizations: finite element ice models, for example.In this case, the concept of exchange grid is not meaningful, since exchange grids are by definition L0.
Instead of the exchange grid, one can choose G to be just about any L0 grid, of resolution at least as high as the ice mesh, where grid cells do not cross atmosphere cell boundaries: we will call this a generalized exchange grid.The user would then need to develop an appropriate conservative transformation for G → I , whereas the transformation for I → G would remain an area-weighted remapping (see Appendix C).Details of this scheme for a particular ice mesh parameterization are left to the reader.
As in Sect.9, the transformation G → I would cause some numerical dispersion, due to the geometry of the ice mesh.Additional dispersion would be introduced because of mismatch between the L0 grid G and the higher order mesh I -although this could be reduced by making G finer.The trade-offs of choosing G to be a generalized exchange grid vs. the ice grid are the same as in Sect.9: numerical dispersion is introduced in the ice coupling time step, in exchange for simplifying implementation and eliminating dispersion in the atmosphere time step.

Regridding examples: local RM
Having shown how to compute all five required regridding transformations, we now demonstrate them working within a realistic GCM context.We set up a test using the GISS 2 • × 2 1 2 • atmosphere grid (Schmidt et al., 2006), overlapping with the SeaRISE 5 km L0 grid (Bindschadler et al., 2013) and 40 elevation points, spaced every 100 m from 0 m to 4000 m.We chose the exchange grid as G (see Appendix A and Fig. 3 for notation conventions).We ran GISS ModelE in this configuration, using fixed sea surface temperatures corresponding to the years 1996-2005.The ice surface model was run on the elevation grid, producing monthly averages of SMB over Greenland, which we label f E (Fig. 14a).
As expected, the broad pattern shows strong melting in narrow bands at low elevations, along with weak accumulation at high elevations.Atmosphere grid cell boundaries are prominent because precipitation -the primary source of positive SMB -is not downscaled to sub-grid resolution (Leung and Ghan, 1998).The small oscillations at high elevation are artifacts introduced by ModelE's ice surface model.

Example: elevation to atmosphere grid
Figure 14b shows the original field f E computed on the elevation grid and regridded to the atmosphere grid A, via the transformation f A = RMf E .Note how low-elevation ablation regions become broader and weaker (compared to Fig. 14a), whereas the interior of the ice sheet remains about the same.This demonstrates the benefits obtained through the use of elevation points, as compared with running the ice surface model on the atmosphere grid (van den Broeke et al., 2008).
By definition, the transformation RM is conservative on A: we can evaluate conservation properties of other transformations by comparing to Fig. 14b.

Example: elevation to ice grid
Figure 14c shows the original field f E regridded to I , via f I = XMf E using z interpolation for M.This plot looks like a smoothed version of Fig. 14a.Atmosphere grid cell boundaries are still visible because z interpolation is local.
The transformation X : G → I (Fig. 13) introduces a small amount of nonlocality.And although it is conservative over the ice sheet in general, it is not conservative over A. This can be seen (Fig. 15) by regridding f I to A and comparing with the results of Fig. 14b.This plot quantifies the amount of numerical dispersion the simulation will encounter every month when preparing SMB input for the ice model.This dispersion is caused by ice grid cells that overlap more than one atmosphere grid cell.
In most areas, numerical dispersion is low, less than 1 %.That is because most atmosphere grid cells are overlapped by many ice grid cells, with only a few lying on an atmosphere grid cell boundary.However, numerical dispersion can be significant for atmosphere grid cells that just nick the edge of the ice sheet, where a high proportion of their overlapping ice grid cells also overlap with a neighboring atmosphere grid cell.
The numerical dispersion encountered here is inherent in the regridding problem itself, rather than our approach to that problem.Short of using an ice grid whose grid cells do not overlap atmosphere grid cells, we see no obvious way to eliminate this issue.However, we do not believe it to be a serious problem: the total area of ice sheet affected by it will be small.It is important to point out that although X is not (quite) conservative over A, it is still conservative over the ice sheet in general.We have verified this numerically in our examples.

Example: ice to elevation grid
Since we have not yet developed the surface boundary condition described in Sect.2.1, we do not have realistic fields on I to try regridding to E. We will therefore test I → E at this point using f I = XMf E as input.We test the transformation in two steps.First, we test G → E using f G = Mf E as input.Then, we test I → E using f I = XMf E as input.This allows us to evaluate the G → E transformation separately from confounding dispersion factors involved in G → I .
We computed Theoretically, f E and f E should be equal.The conservation property imposed by the quadratic program held: we found that f E ≡ A f G ≡ A f E to machine precision.However, we found measurable differences between f E and f E (Fig. 16).Although differences are small in most cases, one area has a difference of more than 8 out of 37 mm day −1 , resulting in no more than one digit of precision.Differences tend to be 65 Fig. 16.Example of numerical dispersion caused by an ill-posed quadratic optimization problem.July SMB on elevation grid E was regridded to the exchange grid.The reverse transformation was then used to recover the original SMB on E. Plot shows the difference between the result and the original.In theory, the two should be exactly the same.Although mass is conserved, differences appear because the associated quadratic optimization problem is underdetermined in some areas.
greatest in sparsely populated atmosphere grid cells on the edge of the ice sheet.We conclude that the answer to the quadratic program posed in Sect.8.5 is poorly constrained.A different numerical solver than the one we used might in theory result in a better match between f E and f E .
But this is not a problem in practice: the goal of G → E is to obtain a physically plausible field on E with the correct conservation properties.The preliminary tests presented in this section are consistent with that goal.
The transformation I → E is constructed by transforming I → G → E (Sect.9).We used this to compute f E = [I → E](f I ) where f I = XMf E .The conservation property imposed by the quadratic program held: we found that f E was equivalent to f I on A to machine precision.
However, differences between f E and f E are even greater in this case: compare Fig. 17  xample of numerical dispersion in the transformation I → E, caused by ice grid cells that ore than one atmosphere grid cell.July SMB on elevation grid E was regridded to the ice grid erse transformation was then used to recover the original SMB on E. Plot shows the difference he result and the original.Although mass is conserved, differences are due mainly to ice grid overlap more than one atmosphere grid cell, making it impossible to recover the exact original same color scale is used as in Fig. 16. 67

Fig. 17. Example of numerical dispersion in the transformation
I → E, caused by ice grid cells that overlap more than one atmosphere grid cell.July SMB on elevation grid E was regridded to the ice grid I .The reverse transformation was then used to recover the original SMB on E. Plot shows the difference between the result and the original.Although mass is conserved, differences are due mainly to ice grid cells that overlap more than one atmosphere grid cell, making it impossible to recover the exact original SMB.The same color scale is used as in Fig. 16.
original f E .Differences are greatest for elevation points near an atmosphere grid cell boundary.

Regridding examples: almost local RM matrix
In Sects.8.3 and 9, we discussed the trade-offs between using the ice vs. exchange grid for the interpolation grid G. Having shown in Sect. 10 an example of our transformations using the exchange grid as G, we now show how things change if the ice grid is chosen for G by re-running the above examples.In this case, the RM matrix will be mostly local, except for nonlocality caused by ice grid cells that intersect more than one atmosphere grid cell.We say that RM is almost local.
In general, differences in results here vs. those in Sect. 10 are insignificant.However, the choice of I as interpolation grid does change the basis functions used for E, yielding a system in which the transformation E → I is now conservative over A.
With a change in E → I , the transformation E → A is also changed.We calculated the difference between f A computed using this method vs. f A computed in Fig. 14b: this difference is exactly the same as Fig. 15.In practice, this difference does not present a real problem: it is impossible to say which version of f A is more "correct."Both maintain conservation over the ice sheet.As predicted in Sect.9, gridrelated dispersion in I → E is eliminated.Thus, the gridrelated "errors" demonstrated in Fig. 17 are eliminated.
Finally, the use of a nonlocal RM matrix requires use of the algorithm described in Sect.8.1 for A → E. Fig. 18a shows a July precipitation field p A produced by ModelE, for the same month as the SMB fields above.ModelE currently does not downscale precipitation to sub-grid resolution (Leung and Ghan, 1998): it therefore yields a precipitation field on A.
We then computed p E = [E → A](p A ), which could serve as input to the land surface model.Recall that p A is the most physically self-consistent value for p E , but that we choose a different value for p E in order to maintain conservation.
Figure .18b shows the difference between the two, p E − p A .Differences are typically in the range [−0.1, 0.1] (about 2 %).They tend to be relatively constant within each atmosphere grid cell and show no discernible pattern between grid cells.In particular, differences are not related to elevation.Any errors introduced by this scheme will be dwarfed by other precipitation errors in the model.
In this example, we used an L0 ice grid to generate a prototypical almost-local RM matrix.We expect similar results when using an L1 ice grid because the nonlocality in the L0 case was caused by a small number of ice grid cells that overlap more than one atmosphere grid cell.A similar situation exists with any type of higher order mesh.
We conclude by considering, in the case of an L0 ice grid, whether I or G is a better choice for the interpolation grid.In many cases, the use of a nonlocal RM matrix requires significantly more effort in GCM model development.However, our experience shows fewer "surprises" in the transformations when using I as the interpolation grid.In the end, we expect either choice to yield serviceable results that conserve mass and energy in the regridding.

Regridding examples: nonlocal RM matrix
In the above sections, we considered how the transformations in this paper work when the RM matrix is local or almost local.Here, we consider the case in which RM is significantly nonlocal -for example, if bilinear interpolation were used for the transformation E → I .We would expect the "corrections" and dispersive properties that appear in the case of an almost-local RM matrix to be significantly larger.We encountered numerous problems when we attempted the use of a significantly nonlocal RM matrix, constructed using bilinear interpolation for E → I .Most significantly, the A → E transformation produced elevation points of negative precipitation when regridding a typical precipitation field over Greenland (Fig. 19) -an artifact that we believe would be unacceptable to the majority of modelers.It might be possible to find a solution to this problems.In the meantime, it is simpler to avoid nonlocal interpolation schemes such as bilinear interpolation.

Independent elevation grid
The elevation grid E is constructed by adding elevation points to each grid cell of an underlying horizontal grid, which we will denote E 0 .So far we have assumed that E 0 = A, meaning the elevation grid E is derived from the atmosphere grid A. Some GCMs use a horizontal layout for E that is not related to A. In this section, we extend our methods to address that case.
The first problem is a choice of the interpolation grid G.The idea of locality in RM no longer makes sense when E 0 = A. For that reason, G = I is the right choice, X and X will both become identity transformations.The transformation E → I does not need to change, nor does I → A: none of these rely on any specific relationship between E 0 and A. Similarly, E → A = RM can be computed as before.
The only other thing that must change is the construction of A → E: the intuitive definition for (Sect.8.1) no longer makes sense.Instead, we construct an intuitive regridding operation F as follows.Given f A , first regrid to f E 0 using area-weighted remapping.Then convert f E 0 to f E using a "repeat" operator akin to above.We can now apply the quadratic programming-based regrid operator developed in Sect.8.2, using F instead of in Eq. ( 15).
The use of E 0 = A introduces numerical dispersion into the system.This can be seen by evaluating the locality of the round-trip transformation E → A → E. This numerical dispersion is a fundamental consequence of the use of a different underlying grid for the ice surface and atmosphere models.

Regridding in elevation space
We have tacitly assumed so far that there is one single fixed elevation grid E with one fixed set of basis functions.This is not the case in a changing climate, since the basis functions used for E depend on ice elevation and extent (Sect.7.1).The user might wish to explicitly move elevation points as well -for example, to track a mountain glacier as it moves up or down.When the vector space E changes, the ice surface Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | mplausible precipiation field produced by A → E when bilinear interpolation is used.July on (Fig. 18a), was regridded to the elevation grid E using the A → E reverse transformation, matrix for E → A is constructed using (non-local) bilinear interpolation.Note the unphysical egative precipitation) for some elevation points!69 Fig. 19.Implausible precipitation field produced by A → E when bilinear interpolation is used.July precipitation (Fig. 18a) was regridded to the elevation grid E using the A → E reverse transformation, where the matrix for E → A is constructed using (nonlocal) bilinear interpolation.Note the unphysical artifacts (negative precipitation) for some elevation points.model state must be regridded from the old elevation grid to the new elevation grid.We address that issue in this section.
Assume two elevation grids, an "old" grid E and a "new" grid F .We wish to conservatively regrid a field f E to f F .In this case, f will not be a flux variable, but rather a conserved state variable of the ice surface model: snow depth, water fraction, etc.
This problem can be approached by examining the system of grids and transformations available in Fig. 13 when one has multiple elevation grids (see Fig. 20).By "following the arrows," the most direct way to transform f E to f F is to first regrid to the interpolation grid, computing f G = M E f E .Then use the procedure in Sect.8.5 to compute f F .Because all these transformations are conservative on the atmosphere grid A, the end result f F will be equivalent to f E on A.

Conserved model state
In previous sections, f represented fluxes between models, which are conserved.In this section, f is a model state variable of the ice surface model.In order for this regridding procedure to be physically meaningful, the model state must be expressed in terms of conserved quantities.
Not all ice surface models are formulated in terms of conserved quantities.In this case, the regridding procedures may still be used, as long as the ice flow model can be converted to/from a form that is expressed in conserved quantities.For example, an ice surface model might track the temperature, mass and water fraction of the top layer of ice.Temperature is not conserved, so this regridding procedure cannot be used directly.However, model state can be converted to enthalpy and mass alone (Aschwanden et al., 2012).These quantities are conserved and can be correctly regridded with conservative transformations.After regridding to the new elevation grid, model state can then be converted back to the original nonconserved parameterization.

When to regrid
Regridding in elevation space is not just required if one changes elevation points, but also any time that the transformation M for E → I changes.Any time that happens, the regridding described in Sect.14 should be used.Relevant events that change M include 1. changes in ice topography: for example, as an ice sheet inflates or deflates due to changing climate; when an ice grid cell changes elevation, the weights in M used to compute it will change; this change will either be  continuous (as in z interpolation) or will jump at certain thresholds (as in elevation class interpolation).
2. changes in ice sheet extent, as an ice sheet grows or melts; for example, if an ice sheet shrinks, then some ice grid cells will no longer participate in the regridding process, and their associated coefficients in M will turn to zero.
3. changes in the elevation or number of elevation points.

Glint2 coupling library
Here we describe Glint2, an open-source implementation of the transformations developed in this paper.Glint2 is a GCM-ice flow model coupling library whose core function is to compute the five transformations required for a coupled GCM-ice flow model system (Fig. 2).These transformations are computed based on a variety of factors: atmosphere and ice grid geometry, ice topography and extent, and elevation levels chosen for the elevation grid (Fig. 21).Glint2 does not just compute these transformations and provide them as a library, it also provides an application programming interface (API) for GCMs to use to couple with ice flow models.As a realization of the mediator design pattern (Gamma et al., 1994), Glint2 is able to shield the GCM from a number of details of the coupling.Although it uses a different codebase, Glint2 has the same purpose as GLINT (Glimmer Interface; Rutt et al., 2009), namely to build a coupling library between GCM and ice models.The GCM programmer who wishes to couple their GCM with an ice model need only do the following: 1. implement an elevation points scheme, and move the ice surface model to the elevation grid.
2. request the RM matrix from Glint2, and then apply it as needed during normal model run; this might be more complex than it sounds, depending on the GCM's system of domain decomposition; storage and application of the RM matrix is simplified if the GCM author knows in advance that it is local.

ask Glint2 to regrid ice surface model inputs from
A → E as needed; this step is only required if RM is nonlocal; otherwise, the transformation is simple enough for the GCM to do on its own.4. accumulate ice surface fluxes on the elevation grid, and pass them to Glint2 every coupling time step.
5. apply fields returned from the ice flow model via Glint2; these fields are returned on the elevation and atmosphere grids as appropriate, eliminating the need for the GCM to regrid them.
Note that the GCM does not need to know anything about the ice grid or ice flow model.Interface code is added to Glint2, not the GCM, for every ice flow model one wishes to support.Over time, we expect the number of supported GCMs and ice flow models to increase, according to researcher demand.This structure is useful because it will give practitioners a way to try different ice models with a GCM fairly easily.

Adoption issues
We expect that Glint2 could be useful to anyone with a GCM who wishes to couple it with an ice flow model.However, many GCMs have centralized regridding strategies that Glint2 does not really fit into.We do not believe this should be a significant barrier to adoption for two reasons: -In most cases, the elevation grid will have unusual "customized" basis functions (Fig. 12).Centralized GCM regridding schemes are not typically equipped to regrid to/from the elevation grid.Nor are they equipped with the algorithms required for the reverse transformations.
-Glint2 hides all details of the ice grid from the GCM, communicating with the GCM with fields on the atmosphere and elevation grids.Because all ice grid-related issues are encapsulated in Glint2, it should not matter to the GCM how or even whether regridding to/from the ice grid is accomplished.As far as the GCM is concerned, the ice flow model might as well be running on the elevation grid.
Another barrier to adoption is the fact that Glint2 is packaged as a library.Many GCM projects are reluctant to add additional library dependencies, due to the complications such dependencies introduce in the build process.We do not believe this should be a serious problem because coupling with an ice model already involves the use of outside libraries.With or without Glint2, the GCM must still manage additional external dependencies.

Availability
Glint2 is written C++ and designed to couple with GCMs and ice flow models written in Fortran 90/2003, C or C++.It also comes with a Python interface, making it easy to test and plot sample regridding problems before incorporating a coupling strategy into a GCM.Glint2 source and documentation are available for download (Fischer, 2013).

Discussion and conclusions
This paper focuses on a system of conservative regridding strategies needed to support tight two-way coupling between a GCM and an ice flow model in the context of elevationbased downscaling in the GCM.Past efforts at one-way coupling have produced downscaling methods that provide SMB fields that match remarkably well with observations and with regional models.However, these efforts used an inconsistent set of transformations, which would result in nonconservation of mass and energy in a two-way coupled system.
In order to achieve consistency, we began by recognizing the elevation grid as an integral part of the coupling problem, along with the atmosphere and ice grids that had previously been considered.We observed that the ice surface model runs on the elevation grid and that the downscaling step (from elevation to ice grid) is actually another form of regridding.We have therefore transformed a coupling problem involving two models and two grids into one with three models and three grids.
We analyzed the regridding transformations needed in a typical two-way coupling and determined that five such transformations are required: three "forward" transformations and two "reverse" transformations.We then set out to develop a consistent implementation for these five transformations, starting with the forward transformations.
We observe that conservation of mass and energy requires consistency between the transformations from the elevation to ice grid (E → I ) and elevation to atmosphere grid (E → A).We achieve consistency by allowing the user to choose E → I and then constructing a transformation for E → A that is consistent with the user's choice.This is a good approach because it allows the user freedom in choosing a downscaling transformation for E → I .Our transformations imply a set of basis functions for the elevation grid, which we demonstrated in plots.
We then addressed the reverse transformations, using the notions of consistency developed in the previous sections.Problems of underdeterminism and overdeterminism, caused by mismatches in dimensionality between grids, are addressed by using quadratic optimization for these transformations.
The result is a set of five conservative transformations needed to support a two-way coupling of GCMs and ice models.We defined a property of our E → A transformation called locality.We showed a number of theoretical and practical benefits if that transformation is local.We implemented all five transformations and demonstrated that they work in practice on realistic input fields.
Although three of those five transformations are wellknown in the literature (Ramshaw, 1985;Lipscomb et al., 2013), the reverse transformations are new.Most importantly, we showed how to choose a set of grids, basis functions and regridding schemes so that all five regridding operations are conservative.Note that the elevation grid uses a "custom" set of basis functions.This is a significant step beyond past efforts that have focused on conservative regridding to L0 pre-chosen grids (Ramshaw, 1985).
Our downscaling transformation from the elevation to ice grids is taken from previous one-way studies (Lipscomb et al., 2013).Our results "look" almost identical, except for subtle differences in the other related transformations needed to ensure conservation.We therefore provide practitioners with a way to bring already proven downscaling techniques into a conservative two-way coupled setting.
We went on to package these transformations in Glint2, a coupling library.As a piece of software, Glint2 serves as a mediator (Gamma et al., 1994) between the GCM and ice flow model, insulating each from the specific details of the other.This will simplify the coupling of GCMs with multiple ice models, as needed to support future research.
Using the techniques from Appendix D, one can compute the polygon β j as the intersection between element j and region B. Green's Theorem may then be used to compute β j N ij (x, y)dA. (C3) We can then apply our definition of sub-basis functions to achieve our goal: B N i (x, y)dA = j β j N ij (x, y)dA. (C4) This section has provided just an outline of the process.The algebra can become complex at times, and a symbolic computation system such as Maxima can be useful.But in the end, integration of a vector f I over an area B is computed as a linear combination of the elements of f I , just as with L0 grids.

Appendix D Computing the exchange grid
Appendices B and C can be applied repeatedly to compute conservative regridding matrices from an L0 or L1 grid to an L0 destination grid.Both procedures assume a way to compute polygonal intersections between two sets of polygonsalso known as an exchange grid (Balaji et al., 2006).
This task is simple in principle, using modern computational geometry packages such as CGAL (2013).If A and G are two sets of polygons, we explicitly construct each polygon in A and G and then compute pairwise intersections between the two sets (Chin and Wang, 1983).This algorithm provides not just integration formulas, but also the actual polygonal outlines of the exchange grid.
If an appropriately robust polygon intersection algorithm is used, our procedure can deal with nonconvex polygons and other possible irregularities.This is not just a theoretical issue: latitude-longitude grid cells commonly used in GCMs are not convex in spherical or planar geometry.In other cases, practitioners might wish to use grid cells consisting of multiple disjoint polygons.
With issues of polygonal intersection taken care of by prepackaged algorithms, the main challenge here is to find those intersections in a scalable manner.The naive algorithm is to write a nested loop, requiring |A| × |G| iterations.This algorithm, with O(n 2 ) complexity, takes too long even on grids commonly used by GCMs and ice flow models today.Most of the "intersections" of A i and G j will result in nothing, because A i and G j are far from each other and obviously do not intersect.
This problem is solved by using an R-tree (Guttman, 1984) to avoid having to consider intersections of grid cells that are far away from each other.The procedure works as follows: first load all the grid cell outlines of A into the R-tree, indexed by their bounding rectangles.Then loop through each grid cell in G, checking the R-tree for any grid cells in A that it might intersect with.The polygon intersection algorithm is run on each of those grid cell pairs to determine the exact outlines of the exchange grid cells.Running time is cut down to a more reasonable O(n log n).
Past algorithms exist to compute regridding matrices by integrating functions around each cell in an exchange grid (Ramshaw, 1985;Dukowicz and Kodis, 1987;Jones, 1999).These algorithms were originally presented in terms of quadrilateral meshes, but they can also be applied to arbitrary meshes.However, they never explicitly compute the exchange grid polygonal outlines.By explicitly computing an exchange grid and then using that to produce the regridding matrix, we have presented here a procedure that is conceptually simpler, possibly more flexible, but almost certainly not faster to run.

Fig. 1 .
Fig. 1.Configuration of the three models, and the separate spatial domains they occupy.The atmosphere model (A) is coupled to the ice surface model (E), which tracks the top few meters of ice.The ice surface model (E) is then coupled to the dynamic ice flow model (I ).The dynamic ice flow model, operating at long time steps, is insulated from high-frequency surface processes in the iceatmosphere interaction.

Fig. 2 .
Fig.2.Data flow (blue arrows) for the coupling between atmosphere, ice surface and dynamic ice flow models (boxes).Since the three models run on different grids, regridding operations (ovals) are required at each step.Figure21shows the inputs required to compute these regridding operations.

AFig. 3 .Fig. 3 .
Fig. 3. Definition of symbols used throughout the paper.See Appendix A for notational convention

Fig. 4 .
Fig. 4. Interpolated surface mass balance (SMB) function within one atmosphere grid cell.These values come from one month (July) of a run of GISS ModelE with elevation points.The vertical dashed line indicates the mean elevation of the grid cell, "seen" by the atmosphere.Dots and squares represent the results of extrapolated SMB computations at other elevations: dots represent elevations found within the grid cell, whereas squares represent elevations outside the cell's range.SMB values at intermediate elevations are interpolated.

Fig. 7 .Fig. 7 .
Fig. 7. Results of (a) Z Interpolation and (b) bilinear interpolation, to develop a downscaled field on the ice grid from an elevation point field.Bilinear interpolation eliminates the discontinuities present at atmosphere grid cell boundaries.However, total SMB is changed, particularly for localized events such as snowstorms: for example in the Northwest.This could introduce biases in the long-term evolution of the ice sheet, even though both schemes can be made to conserve mass and energy.This figure is for demonstration purposes only.See Sects. 10 and 11 for thorough regridding examples.56 Fig. 7. Results of (a) z interpolation and (b) bilinear interpolation,to develop a downscaled field on the ice grid from an elevation point field.Bilinear interpolation eliminates the discontinuities present at atmosphere grid cell boundaries.However, total SMB is changed, particularly for localized events such as snowstorms, for example in the northwest.This could introduce biases in the long-term evolution of the ice sheet, even though both schemes can be made to conserve mass and energy.This figure is for demonstration purposes only.See Sects. 10 and 11 for thorough regridding examples.

Fig. 8 .
Fig. 8. Setup for bilinear interpolation.The value in an ice grid cell (square) will be the sum of the values at the four nearest atmosphere grid cell centers, weighted by -longitude and -latitude along the axes.

Fig. 9 .
Fig. 9. Traditional elevation class schemes are equivalent to running the ice surface model on an L0 grid, where grid cell outlines are created by the intersection of the atmosphere grid and elevation contours.One such grid is shown in this figure.Note that grid cells extend only as far as the ice sheet; the grey line shows the Greenland coast.

Fig. 10 .
Fig. 10.Interpolated SMB function within one atmosphere grid cell when using elevation classes.The resulting piecewise constant interpolation is almost never preferable to the piecewise linear interpolation in Fig. 4. Traditional elevation class schemes are discouraged because they offer no benefit over first-order z interpolation.

Fig. 12 .
Fig. 12. Unitless basis functions for the elevation grid E, constructed using 20 elevation points and z interpolation (the exchange grid was used as the interpolation grid).The grey box represents one atmosphere grid cell on the west coast of Greenland, with the coastline shown as black lines.White lines in the atmosphere grid cell represent elevation contours corresponding to each elevation point.The basis functions corresponding to elevation points at 950 m, 1150 m and 1350 m are shown.Note that basis functions overlap and are not orthogonal.Because of the z interpolation, each basis function has maximum value at its corresponding elevation, but it has a non-zero support up to one elevation point away.

Fig. 13 .
Fig.13.The five grids used in the coupling problem, and transformations between them.The (linear) interpolation step M, from E to G, is chosen by the user.The (linear) transformation R, from G to A, is an area-weighted remapping step.These two fully constrain the transformation from E to A as RM.If I is L0-parameterized, then G is the exchange grid between A and I , and X and X are areaweighted remapping transformations.If I is L1-parameterized, then G is equivalent to I , making X and X the identity.P is a diagonal rescaling operation between A and A. Reverse transformations required by the coupled system are shown as dotted lines.
Suppose we have a flux field f G on the interpolation grid G.

Fig. 14 .
Fig. 14.July SMB computed by ModelE on the elevation grid E (a), and regridded to the atmosphere grid A via the exchange grid (b) or to the ice grid I using z interpolation (c).For plotting purposes, each elevation point in (a) is assigned to a nearby region of similar elevation; the value of the elevation point is then plotted in its corresponding region.Bold numbers on the color scale indicate the extreme values of the plot.Elevation contours are plotted.

Fig. 15 .
Fig. 15.Difference in July SMB obtained with E → A, using two choices for the interpolation grid G: either the exchange grid or the ice grid.Differences (exchange minus ice) are due to ice grid cells that overlap more than one atmosphere grid cell.
to Fig. 16.This is because the transformation X : G → I is not conservative on A. The end result of transforming E → G → I → G → E will produce an f E that is significantly different from the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 18 .Fig. 18 .
Fig. 18.Non-physical changes in precipiation field introduced by A → E in order to ensure conservation.Panel (a) shows a July precipitation field computed by ModelE on the atmosphere grid.This field was regridded to the elevation grid E using an almost-local RM matrix for E → A. It produced a slightly different result than the physicall intuitive procedure of "repeating" precipitation values for elevation points in each atmosphere grid cell -which is not conservative.Those differences are plotted in panel (b).68

Fig. 20 .
Fig. 20.When ice extent, ice topography or elevation points change, the basis functions for the elevation grid E change along with it.Ice surface model state, which exists on E, must be regridded to the new set of basis functions.Shown is a grid system that can serve as a map for this regridding: E is the old elevation grid, F is the new one and G is the interpolation grid (same for old and new).Ice surface model state may be regridded from E to F by first regridding E → G, then G → F .Note that this diagram is a simplified version of Fig. 13 in which two different elevation grids have been accounted for.

Fig. 21 .
Fig.21.The Glint2 workflow used to compute the regridding operations required by the fully coupled system (Fig.2).Glint2 produces regridding operations based on a variety of factors: atmosphere and ice grid geometry, ice topography and extent, and levels chosen for the elevation grid.Glint2 must recompute the operations when any of these factors changes.