A note on the influence of atmospheric model resolution in coupled climate – icesheet simulations

Coupled climate–ice-sheet simulations have been growing in popularity in recent years. Experiments of this type are however challenging as ice sheets evolve over multi-millennial time scales, which is beyond the practical integration limit for most Earth-system models. A common method to increase model throughput is to trade resolution for computational efficiency (compromises accuracy for speed). Here, we analyze how the resolution of an atmospheric general circulation model (AGCM) influences the simulation quality of a standalone ice-sheet model. Four identical AGCM simulations of the Last 5 Glacial Maximum (LGM) were run at different horizontal resolutions: T85 (1.4◦), T42 (2.8◦), T31 (3.8◦), and T21 (5.6◦). These simulations were subsequently used as forcing of an ice-sheet model. While the T85 climate forcing reproduces the LGM ice sheets to a high accuracy, the intermediate resolution cases (T42 and T31) fail to build the Eurasian Ice Sheet. The T21 case fails in both Eurasia and North America. Sensitivity experiments using different surface mass balance parameterizations improve the simulations of the Eurasian ice-sheet in the T42 case, but the compromise is a substantial ice buildup in Siberia. The 10 T31 and T21 cases are not improving in the same way in Eurasia, though the latter simulates the continent-wide Laurentide Ice Sheet in North America. The difficulty to reproduce the LGM ice sheets in the T21 case is in broad agreement with previous studies using low-resolution atmospheric models, and is caused by a deterioration of the atmospheric climate between the T31 and T21 resolutions. It is speculated that this deficiency may demonstrate a fundamental problem using low-resolution atmospheric models in these types of experiments. 15


Introduction
Experiments with coupled climate-ice-sheet models have become increasingly popular in recent years, much thanks to coordinated international modeling initiatives such as the "Ice Sheet Model Intercomparison Project" (ISMIP6) (Nowicki et al., 2016) and the "Pliocene Ice Sheet Modelling Intercomparison Project" (PLISMIP) (Dolan et al., 2012).This is a particularly challenging enterprise as ice sheets have a high thermal inertia that makes their response time considerably greater than all other components of the climate system-the time scale depends on the application but it typically ranges from 10 3 to 10 5 years.
Simulations of this length are beyond the computational capacity of most Earth-system models, and different techniques to bypassing this problem have therefore been devised.Some of the more popular ones used for simulating the ice-sheet evolution over glacial time scales include: 1 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-235Manuscript under review for journal The Cryosphere Discussion started: 7 November 2017 c Author(s) 2017.CC BY 4.0 License.
(i) Force a standalone ice-sheet model with a transient climate record obtained by interpolating between the climate extremes over the period of interest (often simulations of the pre-industrial and the Last Glacial Maximum; PI and LGM, respectively).The interpolation weights are typically derived from oxygen isotope ratios in Greenland and Antarctic ice cores (e.g.Charbit et al., 2007;Fyke et al., 2014).
(ii) Use an asynchronous coupling between an ice-sheet model and a general circulation model (GCM).The ice-sheet model, which is computationally cheaper than the GCM, is then run multiple years between each update of the model climate (e.g.Liakka et al., 2011;Liakka, 2012;Herrington and Poulsen, 2012;Löfverström et al., 2015).
Although no attempt is made here to assess how these different approaches compare to one another, we conclude that they all rely on a number of assumptions and simplifications that potentially can influence the results.For example: (i) assumes that the glacial climate changes as a linear combination of the PI and LGM states, which is at odds with both modeling and proxy evidence of highly nonlinear circulation changes throughout this time period (e.g., Jackson, 2000;Zhang et al., 2014;Lora et al., 2016;Löfverström et al., 2016Löfverström et al., , 2014;;Löfverström and Lora, 2017); (ii) the asynchronous coupling introduces abrupt changes in the boundary conditions, which may force the climate into an unphysical state at the beginning of each coupling interval; (iii) simplified models often rely on statistical dynamics/physics where almost all interactions are prescribed or represented by first-order linear assumptions.
In addition, one issue that has received little attention in the literature is what role the atmospheric grid resolution-the horizontal mesh on which the model equations are discretized-plays in coupled climate-ice-sheet experiments.Simplified models often utilize a coarse horizontal grid for computational efficiency.For example, the atmospheric component of CLIMBER-2 has a horizontal resolution of approximately 10 • × 51 • (Petoukhov et al., 2000), LOVECLIM runs on a 5.6 • × 5.6 • resolution grid (Goosse et al., 2010), and FAMOUS on a 5 • × 7.5 • grid (Smith et al., 2008).These are to be compared with the 1 • × 1 • operational resolution of many modern GCMs (e.g.Flato et al., 2013).
Although a higher resolution is not automatically synonymous with a better model, it generally means that smaller scale phenomena can be resolved, which in turn reduces the need for explicit (parameterized) diffusion.Note that horizontal diffusion is not only influencing (damping) horizontal motions but it may also impact vertical transport and, e.g., convection (Polvani et al., 2004).Several studies have shown that the numerical convergence breaks down somewhere between the T31 (3.8 • ) and T21 (5.6 • ) resolutions in an atmospheric GCM, which (presumably in part due to an increased diffusion rate; Magnusdottir and Haynes, 1999) degrades the representation of even the largest scale atmospheric phenomena, such as planetary waves and jet streams (Polvani et al., 2004;Magnusdottir and Haynes, 1999;Dong and Valdes, 2000;Löfverström et al., 2016).Note that this resolution limit appears to be independent of model physics; e.g.Polvani et al. (2004) and Löfverström et al. (2016)  Here we investigate how the atmospheric resolution influences the equilibrium state in a standalone ice-sheet model simulation.The climate forcing is derived from four identical LGM simulations run at progressively coarser horizontal grids: T85, T42, T31, and T21 (Table 1).While the T85 climate forcing manages to reproduce the LGM reconstruction to a high accuracy, the intermediate resolution cases (T42 and T31) fail to reproduce the Eurasian Ice Sheet, and the T21 case fails in both continents.These results suggest that a "sufficiently" high atmospheric resolution may be required to ensure the quality of (coupled) climate-ice-sheet model experiments.
The models and experiment design are presented in Section 2, the results from the atmospheric model and the ice-sheet model are described in Sections 3 and 4, followed by a more general discussion in Section 5.

Ice sheet model
We use the three-dimensional ice-sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets, version 3.1), run at a 80 km resolution grid that covers most of the Northern Hemisphere.The model treats ice as an incompressible, viscous, and heat-conducting fluid (Greve, 1997), using the shallow-ice approximation (Hutter, 1983) subjected to Glen's flow law (with stress exponent n = 3) (e.g.Van der Veen, 2013).A Weertman-type sliding scheme is also applied (Weertman, 1964).
We run the model in the so-called "cold-ice mode", which means that temperatures exceeding the pressure melting point are artificially reset to the pressure melting temperature.The global sea level is lowered by 120 m to reflect LGM conditions and marine ice is allowed to form where the bathymetry is less than 500 m, otherwise instantaneous calving is applied.The geothermal heat flux is set to a constant values of 55 mW m −2 over the whole domain, and the bedrock relaxes toward isostatic equilibrium with a timescale of 3 kyrs, assuming a local lithosphere and relaxing asthenosphere (Greve and Blatter, 2009).All simulations started from ice-free conditions (interpolation of atmospheric fields is described in Section 2.3) and were run for 150,000 years to ensure that the ice sheets reach their steady state extent.

Ablation parameterizations
SICOPOLIS uses the positive-degree-day (PDD) method to parameterize ablation.The annual melt-potential is estimated from the integrated sum of positive temperatures each year (Braithwaite and Olesen, 1989;Reeh, 1991), assuming that the daily temperatures are normally distributed about the monthly mean value (Calov and Greve, 2005).
Following Charbit et al. (2013), we test the sensitivity to the surface-mass balance scheme using three different PDD-based ablation models: the default parameterizations in SICOPOLIS (based on Reeh, 1991), plus the ones presented in Fausto et al. (2009), andTarasov andPeltier (2002) (henceforth referred to as SICOdef, FST09, and TP02, respectively).These parameterizations use different methods for calculating the degree-day factors for snow and ice (β snow and β ice , respectively), refreezing fraction of melt water (P max ), and standard deviation (day-to-day variability) of temperature (σ).All these parameters are set to numerical constants in SICOdef (β snow = 3 mm day −1 K −1 , β ice = 12 mm day −1 K −1 , P max = 0.6 and σ = 5 • C), while they take on slightly more elaborate expressions in the other parameterizations (see below).

The FST09 model
The standard deviation of daily temperature (σ) is here assumed to change with elevation at a rate of 1.2224 A similar elevation dependence is also applied to P max to account for the increasing probability of melt water refreezing at higher elevation.No refreezing of melt water (P max = 0) is assumed below 800 m, and total refreezing (P max = 1) above 2000 m.
In addition, the FST09 model uses a temperature dependent degree-day factor for ice that varies from β ice = 7 mm day −1 K −1 for warm summer (June-July-August) conditions (≥ 10 • C), to β ice = 15 mm day −1 K −1 for cold summer temperatures (≤ −1 • C).A cubic change is applied for intermediate temperatures.The degree-day factor for snow is a constant with the same numerical value as in SICOdef (β snow = 3 mm day −1 K −1 ).

The TP02 model
The TP02 model uses a similar temperature-dependent parameterization of β ice as in FST09, but with bounds: β ice = 8.3 mm day −1 K −1 , and β ice = 17.22 mm day −1 K −1 for warm (≥ 10 • C) and cold (≤ −1 • C) summer temperatures, respectively.A parameterization of this kind is also applied to β snow that varies between β snow = 4.3 mm day −1 K −1 and β snow = 2.65 mm day −1 K −1 , respectively.The standard deviation of temperature is set to a constant value of σ = 5.2 • C. The refreezing scheme is also more comprehensive (based on Pfeffer et al., 1991;Janssens and Huybrechts, 2000), including both thermodynamics (latent heat release due to refreezing) and pore trapping components.

Climate evolution in SICOPOLIS
The surface temperature and precipitation (over the evolving ice sheets) are calculated using the method described in Liakka et al. (2016), which in turn is based on the general methodology outlined in Charbit et al. (2002Charbit et al. ( , 2007)).While the temperature decreases linearly with height at a fixed lapse rate γ (= −6.5 × 10 −3 K m −1 ), the precipitation changes exponentially as a function of temperature; see Eqs. 1 and 2 in Liakka et al. (2016).The distribution of liquid and solid precipitation is also assumed to vary with temperature: 100% solid precipitation falls if the monthly mean surface air temperature is below −10 • C, and 100% liquid if it is higher than 7 • C. The distribution changes linearly for intermediate temperatures (Marsiat, 1994).The baseline temperature and precipitation fields are (bilinearly) interpolated from the atmospheric (LGM) climatology, using the above lapse rate to correct for elevation biases due to different grids and horizontal resolutions.

Atmospheric model
The atmospheric climate forcing is produced with the Community Atmospheric Model version 3 (CAM3) (Collins et al., 2004(Collins et al., , 2006b)), using four different spectral (horizontal) resolutions: T85, T42, T31, and T21, corresponding to an approximate grid The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-235Manuscript under review for journal The Cryosphere Discussion started: 7 November 2017 c Author(s) 2017.CC BY 4.0 License.spacing of 1.4 • , 2.8 • , 3.8 • , and 5.6 • , respectively (Table 1).The model uses identical parameterizations (same equations) at all resolutions (Collins et al., 2004(Collins et al., , 2006b), however the model climate is tuned by varying twelve parameters governing the representation of clouds and precipitation (convective and stratiform), biharmonic diffusion, and integration time step in order to satisfy the Courant-Friedrichs-Lewy (CFL) condition; some parameter settings are outlined in Table .1, for a complete description of the model see Collins et al. (2004).Note that the model physics is represented in grid space, while the dynamics is discretized in spectral space.The effective diffusion rate is scale dependent and modulated by (horizontal) wave number in the vorticity and divergence equations (Section 3.1.14in Collins et al., 2004).
The planetary boundary conditions are set to reflect LGM conditions, including the orbital parameters and greenhouse gas concentrations outlined by the Paleoclimate Modeling Intercomparision Project (PMIP) (e.g.Kageyama et al., 2017), the ice sheet reconstruction presented by Kleman et al. (2013) (raised to the height of the ICE-5G reconstruction (Peltier, 2004) to amplify the ice formation in those areas in SICOPOLIS), and prescribed monthly varying sea-surface conditions (sea-surface temperature and sea-ice extent) from the (LGM) simulation in Otto-Bliesner et al. (2006).The vegetation cover in non-glaciated areas is prescribed as the modern distribution.
The grid-resolved boundary conditions were spectrally interpolated from the T85 grid to the coarser resolutions in order to ensure an identical setup.However, intrinsic smoothing in the interpolation process lowers the maximum ice-sheet height by as much as 200 m in the T21 experiment; similar but smaller elevation changes are also found in the T42 and T31 cases.All simulations were run for 12 years, from which monthly climatologies were created over the last 10 years.The relatively small number of years is motivated by the use of prescribed sea-surface conditions.A longer sampling rate may alter details in the climatologies but is not expected to change the first order conclusions from the study.
Figures 1 and 2 show the surface temperature and precipitation climatologies that are used as forcing of the ice-sheet model; the full fields are presented in the left columns (panels a,b,c,d), and the difference with respect to the T85 case are shown in the right columns (panels e,f,g).We focus on the surface temperature in boreal summer (June-July-August; JJA) as this is the primary ablation/melt season, but the cumulative sum of precipitation over the year (total annual) as ice can form in all seasons in regions with a positive surface mass balance.
The JJA surface temperature is to first order similar in the two intermediate resolution cases (T42 and T31; Figs.1e,f), featuring a localized warming with respect to the T85 simulation over the northwestern parts of the Laurentide Ice Sheet, the interior of the Greenland Ice Sheet, and the southwestern parts of the Eurasian Ice Sheet.This is thought to be a response to the lowering (smoothing) of the resolved topography on the coarser grids.There are also regions that experience a net cooling with respect to the T85 case, e.g.along the southern margin of the Laurentide Ice Sheet (in the T31 case), and northeast of both the Laurentide and Eurasian Ice Sheets (in the T42 case).The largest differences in precipitation are found in the tropics (dominated by localized convection), while only small changes are found in the midlatitude storm tracks where large-scale The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-235Manuscript under review for journal The Cryosphere Discussion started: 7 November 2017 c Author(s) 2017.CC BY 4.0 License.precipitation patterns are dominating (Fig. 2).Note that planetary waves are still reasonably well captured at the T31 resolution (Magnusdottir and Haynes, 1999;Löfverström et al., 2016).
The T21 case shows a fairly different response with a considerable warming over most of the world's topography (Fig. 1g), including over the LGM ice sheets.This is likely a response to the lower mean-height of the resolved topography (smoothing from the interpolation process), and possibly also from a general degradation of the model climate (see further discussion in Section 5).The midlatitude precipitation field is also considerable altered with respect to the T85 case, with substantially lower precipitation in the eastern parts of the midlatitude storm tracks and thus over the southwestern parts of the ice sheets (Fig. 2).
Note that this is likely due to the model's inability to resolve planetary waves (and hence individual cyclones) at this resolution (Polvani et al., 2004;Magnusdottir and Haynes, 1999;Löfverström et al., 2016).

Ice sheet model results
The left column in Fig. 3 shows the equilibrium ice-sheet extent when using the default SMB parameterization in SICOPOLIS (SICOdef).The ice sheets forming under the high resolution atmospheric climatology (T85; panel 3a) are in close resemblance with the target extent (indicated by solid contours; Kleman et al., 2013).There is essentially only too much ice extending along the Siberian Arctic coast.
The ice sheets forced by the intermediate resolution climatologies (T42 and T31; panels 3b,c) both adequately reproduce the North American ice sheet, while they fail to build the Eurasian counterpart in agreement with the reconstruction.This one-sided mismatch can be understood from the atmospheric climatologies described in Section 3. The (JJA) warming signal over the southwestern parts of the Eurasian Ice Sheet (Figs. 1e,f) is supposedly the main reason for why ice is not forming in this region.Note that although there is a small reduction of precipitation with respect to the T85 case, the warm surface temperatures are by far the most pronounced feature over the Eurasian Ice Sheet (cf.Figs 1 and 2).These results are in broad agreement with Abe-Ouchi et al. (2013), who showed that the Eurasian Ice Sheet is more sensitive to temperature changes than the North American counterpart.The relatively small temperature change over the Eurasian Ice Sheet is thus sufficiently strong to influence the local ice sheet response there.The warming signal in northwestern North America is located in substantially colder region with a considerably shorter ablation season, and therefore has a smaller influence on the local ice sheet evolution.
The T21 case, on the other hand, struggles to reproduce the LGM ice sheets in both continents.Although ice forms in North America, it fails to build the continent-wide Laurentide Ice Sheet and instead forms two separate ice sheets-a smaller eastern and a larger western dome-separated by a wide gap in the region around Hudson Bay.This response bears some structural similarity to the low-resolution model results shown in Beghin et al. (2014) and Charbit et al. (2013), and also the pre-LGM ice sheets in Calov et al. (2005), and Bonelli et al. (2009).Similar to the T42 and T31 cases, the T21 climate forcing is too warm (and presumably too dry) over the southwestern parts of the Eurasian Ice Sheet area to reproduce the ice sheet reconstruction.
The sensitivity experiments with different SMB parameterizations in SICOPOLIS are presented in Fig. 3e-l.The center column (panels 3e,f,g,h) uses the FST09 ablation model, and the right column (panels 3i,j,k,l) the one described in TP02.
Both these alternative SMB parameterizations help improve the Eurasian ice extent in the T42 case (Figs.3f,j).However, this The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-235Manuscript under review for journal The Cryosphere Discussion started: 7 November 2017 c Author(s) 2017.CC BY 4.0 License.improvement happens at the price of a fairly substantial ice buildup in northern/northeastern Siberia and Alaska, which is areas that were largely ice free at the LGM (Svendsen et al., 2004;Kleman et al., 2013;Löfverström and Liakka, 2016).A broadly similar buildup is also seen in the T85 case when using these SMB parameterizations.
The alternative SMB parameterizations are not improving in the ice-sheet simulations in Eurasia when using the lower resolution climatologies (T31 and T21) (Fig. 3g,h,k,l), but they help the formation of a continent-wide Laurentide Ice Sheet in the T21 case (Fig. 3k,l).

Discussion and conclusions
The results presented in this paper attempt to illustrate, albeit in a highly qualitative way, the importance of a high atmospheric model resolution in coupled climate-ice-sheet experiments.The simplified modeling approach, using prescribed sea-surface conditions from an equilibrated fully coupled LGM simulation (Otto-Bliesner et al., 2006), allows us to effectively isolate the effect of resolution in the atmospheric model.In addition, by prescribing the LGM boundary conditions (sea-surface conditions plus ice sheets) when creating the atmospheric climatologies, we are priming the ice formation to occur in the "correct/desired" areas in the subsequent ice-sheet model experiments.This methodology appears to work well when using the high resolution atmospheric climatology (T85) (see also Liakka et al., 2016), but is less successful when using the climatologies from the lower resolution simulations (T42, T31, and T21; Fig. 3).The difficulty in reproducing the Eurasian Ice Sheet with the lower resolution climatologies suggests that a "sufficiently" high atmospheric resolution is required to ensure the quality of these types of experiments.The precise limit when the atmospheric resolution starts to become an issue is likely model dependent and contingent on a range of different factors, including the SMB parameterization in the ice-sheet model (some evidence of this is shown in Fig. 3).
One piece of information that is rarely mentioned in the paleo-modeling literature is that most Earth-system models are tuned to reproduce the climate of the instrument era (∼1850 to present).Although these models are extremely valuable tools for exploring other time periods as well, it generally means that inter-model discrepancies tend to increase under more extreme forcing scenarios, e.g.glacial conditions (e.g.Braconnot et al., 2007).The results presented here suggest that the model spread could be further exacerbated by a general degradation of the model climate at lower horizontal resolutions.This encapsulates the challenges of coupled climate-ice-sheet simulations (long simulations in general), as many models often trade resolution ("accuracy") for computational efficiency (speed) in order to run transient experiments over glacial timescales.
The atmospheric model used here has been tuned and extensively tested at the T85, T42, and T31 grids (e.g.Collins et al., 2006a;Yeager et al., 2006).However, the T21 resolution only has "functional support", which means that boundary conditions are provided but the model climate has not been tuned to the same standard as the other resolutions (the resolution dependent tuning parameters are broadly the same as in the T42 case).This is probably at least a partial explanation for the apparent degradation of the model climate at T21, though it is possible that this also manifests a more general breakdown of the numerical convergence that appears to happen around this resolution in other models as well (e.g.Polvani et al., 2004;Magnusdottir and Haynes, 1999;Dong and Valdes, 2000).A curious aspect of the T21 experiment is that it fails to build the continent-wide Laurentide Ice Sheet in North America when using the default SMB parameterization (SICOdef).Instead it builds two spatially disconnected ice sheets, with a larger dome (center of mass) on the western side of the continent (Fig. 3d).Several coupled climate-ice-sheet experiments with a low-resolution atmospheric model have shown qualitatively similar results, see e.g.: Calov et al. (2005); Charbit et al. (2013); Beghin et al. (2014).The common denominator for these studies is that they all used CLIMBER-2 to produce the atmospheric forcing fields.We stress that it is not our intention to single out this particular model, but it appears to suffer from similar deficiencies as our T21 case and may therefore help us understand some of these results.In the aforementioned papers the ice sheet tends to be limited to the western/northwestern side of the North American continent (e.g.Charbit et al., 2013;Beghin et al., 2014), little or no ice is established in western Eurasia (e.g.Calov et al., 2005;Charbit et al., 2013;Beghin et al., 2014), and attempts to remedy these shortcomings typically result in substantial ice formation in Siberia and Alaska (see Charbit et al., 2013, who tested the sensitivity to the same PDD-based SMB parameterizations as were used in this study).
These results appear to be largely independent of both the choice of ice-sheet model (the above studies used SICOPOLIS and GRISLI), and the complexity of the SMB parameterization (Charbit et al., 2013;Bauer and Ganopolski, 2017).Although it is not completely fair to compare CLIMBER-2 to a somewhat poorly tuned low resolution version of CAM3 (these models are extremely different in terms of both complexity and general purpose), it is possible that these similarities demonstrate a fundamental problem with low-resolution climate models, which cannot simply be improved by tuning or increasing the model complexity.
found a similar limit using a dry primitive equation model (no model physics) and a comprehensive atmospheric circulation model (fairly sophisticated physics), respectively.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-235Manuscript under review for journal The Cryosphere Discussion started: 7 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .Figure 2 .
Figure 1.Summer (JJA) surface temperature [ • C] from the different resolution atmospheric climatologies.The full fields are shown in the left column (panels a,b,c,d), and the difference with respect to the T85 case is shown in the right column (panels e,f,g).The 500 m ice thickness from the LGM reconstruction is indicated by the heavy contours (interpolated to the different horizontal resolutions).

Figure 3 .
Figure3.Equilibrium ice thickness [m] when using different ablation parameterizations in the surface mass balance scheme: (left) default method in SICOPOLIS; (middle) method byFausto et al. (2009); and (right) method byTarasov and Peltier (2002).The atmospheric climatologies from the (a,e,i) T85; (b,f,j)T42; (c,g,k) T31; and (d,h,l) T21 resolution simulations were used, respectively.The 500 m ice thickness from the LGM reconstruction is indicated by the heavy contours (interpolated to the different horizontal resolutions).

Table 1 .
Selected parameter settings at different spectral resolutions.The corresponding grid resolution [ • ] is outlined in the top row, integration time step [s] in the middle row, and horizontal (biharmonic) diffusion coefficient [10 15 m 4 s −1 ] in the bottom row.