The Glimmer Community Ice Sheet Model (Glimmer-CISM) has been implemented in the Community Earth System Model (CESM). Glimmer-CISM is forced by a surface mass balance (SMB) computed in multiple elevation classes in the CESM land model and downscaled to the ice sheet grid. Ice sheet evolution is governed by the shallow-ice approximation with thermomechanical coupling and basal sliding. This paper describes and evaluates the initial model implementation for the Greenland Ice Sheet (GIS). The ice sheet model was spun up using the SMB from a coupled CESM simulation with preindustrial forcing. The model's sensitivity to three key ice sheet parameters was explored by running an ensemble of 100 GIS simulations to quasi equilibrium and ranking each simulation based on multiple diagnostics. With reasonable parameter choices, the steady-state GIS geometry is broadly consistent with observations. The simulated ice sheet is too thick and extensive, however, in some marginal regions where the SMB is anomalously positive. The top-ranking simulations were continued using surface forcing from CESM simulations of the twentieth century (1850–2005) and twenty-first century (2005–2100, with RCP8.5 forcing). In these simulations the GIS loses mass, with a resulting global-mean sea level rise of 16 mm during 1850–2005 and 60 mm during 2005–2100. This mass loss is caused mainly by increased ablation near the ice sheet margin, offset by reduced ice discharge to the ocean. Projected sea level rise is sensitive to the initial geometry, showing the importance of realistic geometry in the spun-up ice sheet.
Most global climate models have not included interactive ice sheet models, in part because ice sheets were once thought to evolve only on time scales of centuries to millennia. Recent observations, however, have shown that ice sheets can respond to ocean and atmospheric changes on annual-to-decadal time scales and that mass loss from the Greenland and Antarctic Ice Sheets has increased since the 1990s. In Greenland, atmospheric warming has led to greater mass loss through increased surface melting and sublimation (van den Broeke et al. 2009). For both ice sheets, mass discharge through fast-moving outlet glaciers and ice streams has risen (Howat et al. 2007; Rignot et al. 2008), possibly triggered by intrusions of warm seawater beneath floating ice shelves (Shepherd et al. 2004; Holland et al. 2008; Pritchard et al. 2012). During the past decade, the two ice sheets have lost mass at a combined rate of 300 to 500 Gt yr−1 (Rignot et al. 2011), raising the global-mean sea level by an average of ~1 mm yr−1. The recent ice sheet contribution to global-mean sea level rise (SLR) is comparable to the estimated contributions from ocean thermal expansion and melting of smaller glaciers and ice caps (Church et al. 2011) and is likely to increase as the earth's climate warms further.
Estimates of twenty-first-century ice sheet mass loss are highly uncertain. The Intergovernmental Panel on Climate Change Fourth Assessment Report (IPCC AR4) (Meehl et al. 2007) projected 0.18 to 0.59 m of sea level rise by 2100 but excluded ice sheet dynamical feedbacks because existing ice sheet models were deemed inadequate. Semiempirical models (e.g., Rahmstorf 2010; Jevrejeva et al. 2010) typically estimate twenty-first-century sea level rise of ~0.5 to 2.0 m, based on the assumption that the rate of SLR is linearly proportional to changes in global-mean temperature or radiative forcing. The differences between process-based and semiempirical projections have not been fully explained. Improved process-based models, especially of ice sheets and their interactions with the ocean and atmosphere, are needed to better bound the projections.
Ice sheet modelers generally agree on the need for 1) Stokes or higher-order ice flow models with a unified treatment of vertical shear stresses and horizontal plane stresses, 2) grid resolution ~1 km or less near grounding lines (the boundaries between grounded and floating ice) and in other regions where the flow varies rapidly on small scales, and 3) more realistic treatment of physical processes such as basal sliding, subglacial water transport, and iceberg calving (Little et al. 2007; Lipscomb et al. 2009). These improvements are beginning to be incorporated in a new generation of ice sheet models (Gagliardini and Zwinger 2008; Pattyn et al. 2008; Price et al. 2011; Larour et al. 2012; Cornford et al. 2013). The new models are providing improved estimates of ice sheet mass loss and resulting SLR, for example, as part of the Sea-level Response to Ice Sheet Evolution (SeaRISE) and ice2sea multimodel assessment projects (Bindschadler et al. 2013; Shannon et al. 2013).
Improved ice sheet models are necessary, but not sufficient. These models must be coupled to climate models to capture interactions with the atmosphere and ocean. For example, surface ablation may be underestimated if the ice sheet surface does not interact with the atmosphere (Pritchard et al. 2008). At ice sheet margins, floating ice shelves are vulnerable to warming of subshelf ocean waters (Pritchard et al. 2012), but these interactions are only beginning to be modeled (Hellmer et al. 2012). Also, changes in ice sheet topography and surface runoff could alter atmospheric and oceanic circulation patterns (Ridley et al. 2005; Hu et al. 2009). Coupling ice sheet models to climate models is a long-term effort requiring major software engineering changes (e.g., to allow the land–atmosphere and land–ocean boundaries to change in time), fine grid resolution (e.g., to simulate grounding-line migration), and detailed validation studies.
In one of the first ice sheet coupling efforts, Ridley et al. (2005) coupled a 3D thermomechanical ice sheet model to the third climate configuration of the Met Office Unified Model (HadCM3) climate model and used this model to explore long-term deglaciation of the Greenland Ice Sheet (GIS) under elevated CO2. Driesschaert et al. (2007) also studied future melting of the GIS using a coupled model. Subsequently, Ridley et al. (2010) explored reglaciation thresholds and GIS hysteresis, and Gregory et al. (2012) used a coupled ice sheet–climate model to study glacial inception. Each of these models used variations of a positive-degree-day (PDD) model to compute the ice sheet surface mass balance (SMB, the difference between annual accumulation and ablation). In PDD models, surface ablation is computed based on empirical relationships between temperature and melting, rather than detailed surface energy balance calculations. Mikolajewicz et al. (2007) and Vizcaíno et al. (2010) used a more physical energy balance approach to compute the ice sheet SMB in a coupled framework and to simulate evolution of the Greenland and Antarctic Ice Sheets under future warming. Fyke et al. (2011) coupled an ice sheet model to an intermediate-complexity climate model using an energy balance scheme and studied the effects of albedo, refreezing, and climate model bias on GIS geometry in various climate states.
The steady-state ice sheet geometry in a coupled ice sheet–climate model is sensitive to SMB biases (e.g., Pollard 2000; Fyke et al. 2011). For this reason, most coupled models use some form of anomaly coupling to remove first-order climate biases from the surface temperature or SMB. This approach is acceptable for small perturbations from the present state, but may be problematic for larger changes because of the assumption that climate model biases will remain fixed in a changing climate. (A similar argument applies to flux adjustments in coupled atmosphere–ocean models.) Ideally, ice sheet models should be forced with an uncorrected SMB, but generating a surface mass balance that results in realistic ice sheet geometry remains a challenge for large-scale climate models.
Here we describe initial efforts to include dynamic ice sheets in the Community Earth System Model (CESM; http://www.cesm.ucar.edu), a global climate model that is being used extensively for national and international climate assessments (Hurrell et al. 2013). CESM is the successor of the Community Climate System Model (CCSM) (Gent et al. 2011), which consists of four physical models (atmosphere, ocean, land, and sea ice) linked through a coupler. CESM adds a dynamic ice sheet model—the Glimmer Community Ice Sheet Model (Glimmer-CISM)—as a fifth physical component. The combined system, known as CESM1(CISM), also includes a new surface mass balance scheme for ice sheets in the CESM land surface component, the Community Land Model (CLM), along with infrastucture for coupling CLM to Glimmer-CISM (Fig. 1). Work to date has focused on simulating and validating the SMB and steady-state geometry of the Greenland Ice Sheet as a step toward fully interactive ice sheet–climate modeling. This paper describes the model implementation and evaluates the ice sheet response to CESM surface forcing. Two companion papers (Vizcaíno et al. 2013a,b) analyze the GIS surface mass balance in coupled preindustrial twentieth-century (1850–2005) and twenty-first-century (2005–2100) CESM1(CISM) simulations.
The paper is organized as follows. Section 2 describes the ice sheet model, the calculation of ice-sheet surface mass balance in CLM, and the coupling of CLM to Glimmer-CISM. Section 3 describes model spinup, analyzes the steady-state GIS geometry under preindustrial forcing, and uses an ensemble of simulations to identify optimal ice sheet model parameters. Section 4 estimates GIS mass loss and the resulting sea level rise under historical and projected climate forcing. Finally, section 5 presents conclusions and discusses future model development, including more realistic treatment of ice flow in Glimmer-CISM, improved downscaling and snow/ice physics in CLM, and interactive coupling between Glimmer-CISM and other CESM components.
2. Model description
CESM1(CISM) simulates ice sheets using three new model components:
Glimmer-CISM, a dynamic ice sheet model, which computes ice velocities and the resulting evolution of ice sheet geometry and temperature;
a surface mass balance scheme in the Community Land Model, which computes accumulation and ablation at the upper surface of ice sheets; and
software for coupling CLM and Glimmer-CISM.
a. The dynamic ice sheet model
Glimmer-CISM is a thermomechanical ice sheet model that solves three-dimensional equations for conservation of momentum, mass, and internal energy (Rutt et al. 2009). It uses open source code written primarily in Fortran 90. A standalone version of Glimmer-CISM is available from the BerliOS repository (http://developer.berlios.de/projects/glimmer-cism/), and a CESM developmental version is maintained at the National Center for Atmospheric Research. The version initially implemented in CESM is Glimmer-CISM 1.6.
The model's dynamical core, known as Glide, is based on the shallow-ice approximation (SIA). The SIA is valid for slow flow dominated by vertical shear but not for fast flow in ice shelves, ice streams and outlet glaciers, where lateral and longitudinal stresses are important. More sophisticated higher-order models (Pattyn et al. 2008) are needed to simulate fast flow in these regions. The SIA is a reasonable simplification for much of the Greenland Ice Sheet, which has little floating ice, but would be inappropriate for the large ice shelves that border the Antarctic Ice Sheet (AIS). For this reason we have not yet carried out AIS simulations in CESM.
The Glimmer-CISM climate model interface, Glint, was modified to receive an externally computed SMB from CESM. Glint also includes a positive-degree-day scheme that is not used by CESM. Within CESM, Glimmer-CISM supports coupling to a dynamic Greenland Ice Sheet on a 5-km rectangular mesh, based on a polar stereographic projection. (The model can also be run on 10-km and 20-km meshes.) The ice sheet is initialized using ice thickness and topography from the Bamber et al. (2001) dataset, supplemented by high-resolution data around Jakobshavn Isbrae (Gogineni 2012).
The simulations described in this paper use a simple parameterization of basal sliding. In regions where the bed is at the ice melting temperature, the basal ice velocity is linearly proportional to the driving stress, with a constant of proportionality that does not vary in space or time. Where the bed is below the melting temperature, the basal velocity is zero. Calving also is treated simply; the ice thickness is set to zero wherever the ice is found by Archimedes' principle to be afloat. Thus, there are no ice shelves in these simulations.
b. The ice sheet surface mass balance
The ice-sheet surface mass balance is computed by CLM using an energy balance scheme. In such a scheme, the energy available for surface melting is determined from the sum of radiative, turbulent, and conductive heat fluxes. Energy balance schemes are more expensive and complex than positive-degree-day schemes, but also are more realistic. PDD schemes are not ideal for climate change studies because empirical degree-day factors could change in a warming climate. Comparisons of PDD and energy balance schemes (Bougamont et al. 2007) suggest that PDD schemes may be overly sensitive to warming temperatures.
The Community Land Model typically runs on an ~1° grid, which is too coarse to capture topographic variations on the steep slopes of ice sheets. To improve accuracy, the surface mass balance is computed for a specified number of elevation classes in each CLM grid cell that overlaps the Glimmer-CISM ice sheet grid. By default there are 10 elevation classes; the boundaries (in meters) defining these classes are 0, 200, 400, 700, 1000, 1300, 1600, 2000, 2500, 3000, and 10000. Fyke et al. (2011) used a similar scheme to compute the ice sheet SMB in the University of Victoria (UVic) Earth System Climate Model.
There are several advantages to computing the surface mass balance in CLM rather than Glimmer-CISM.
It is computationally cheaper to compute the SMB on the land grid for 10 elevation classes per grid cell than on the 5-km ice sheet grid, which has about 100 times as many grids cells as CLM (at 1° resolution) per unit area of Greenland.
We can use the sophisticated CLM snowpack model (Oleson et al. 2010) instead of maintaining a separate scheme for Glimmer-CISM.
The atmosphere can respond during runtime to ice sheet surface changes, capturing critical albedo feedbacks (Pritchard et al. 2008).
It is easier to ensure that the change in ice mass due to freezing and melting in CLM is equal to the mass change in Glimmer-CISM. (Because of approximations during grid interpolation, mass is not conserved exactly in the current model, but for two-way land–ice sheet coupling it will need to be exact.)
A surface mass balance scheme in CLM can be used not only for ice sheets but also for mountain glaciers and ice caps, whose evolution will be modeled in future versions of CESM.
The CLM hierarchical data structure makes it fairly straightforward to implement multiple elevation classes for a given land surface type, or “landunit.” In the standard version of CLM (without elevation classes), each grid cell is divided into one or more of five landunits: vegetated, lake, wetland, urban, and glacier. Each landunit consists of a specified number of columns. Urban landunits contain five columns, while vegetated, lake, wetland, and glacier landunits contain one column each. Every column in a grid cell receives the same temperature, specific humidity, and downward radiative fluxes from the atmosphere, but each column has its own sensible and latent heat fluxes (computed using Monin–Obukhov similarity theory), upward radiative fluxes, and subsurface temperature and water content. When CLM is configured for ice sheets, glacier landunits are replaced by glacier_mec landunits, which have similar physics but contain a separate column for each elevation class. (Here, “mec” stands for multiple elevation classes.)
The CLM grid cells are initialized with a surface topography consistent with the ice sheet model. That is, the Greenland digital elevation map (Bamber et al. 2001) used to initialize Glimmer-CISM is also used to compute the fractional glacier coverage in each elevation class of each grid cell in the part of the CLM domain that overlaps the ice sheet domain. As the ice sheet evolves, its topography can diverge from the CLM topography, and Glimmer-CISM may require surface forcing at elevations that have zero glacier coverage in CLM. For this reason, CLM does a prognostic calculation for all 10 columns in each Greenland grid cell, not just the columns with nonzero area. Columns with zero area, known as virtual columns, do not affect the atmosphere (since their surface fluxes are given zero weight when summed over a grid cell) but exist only to provide forcing, as needed, for the ice sheet model.
In CLM grid cells with glacier_mec landunits, the atmospheric surface temperature and specific humidity are downscaled from the mean gridcell elevation to the various column elevations using a globally uniform lapse rate (typically 6°C km−1) and holding relative humidity constant. At a given time, the lower-elevation columns in a grid cell can undergo surface melting while columns at higher elevations remain frozen. This results in a more accurate simulation of summer melting, which is a nonlinear function of air temperature. Radiative fluxes and precipitation are not currently downscaled, in part because it is difficult to do so while conserving energy and water. In particular, the partitioning of precipitation between rain and snow is determined by the atmosphere model; rain is not converted to snow (or vice versa) as a function of downscaled temperature because this would violate energy conservation.
CLM has a sophisticated multilayer snow model (Oleson et al. 2010) that allows downward percolation and refreezing of surface meltwater. When meltwater reaches the snow–ice interface or exceeds the holding capacity of the snowpack, it runs off. Runoff is sent to the ocean using the CESM river routing model. As snow ages, it is compacted by destructive metamorphism, overburden pressure, and melting/refreezing. Snow albedo and radiative transfer within each layer are computed using the Snow and Ice Aerosol Radiative Model (SNICAR) (Flanner and Zender 2006). The snow albedo depends on zenith angle, effective grain size, and the presence of aerosols (black carbon, organic carbon, and mineral dust) deposited from the atmosphere. The albedo of bare ice, however, is simply set to a prescribed value in the visible and near-infrared bands. In the simulations described below, the visible and near-IR ice albedos are 0.60 and 0.40, respectively, in approximate agreement with observed values for melting glacier ice (e.g., Bøggild et al. 2010).
In the standard version of CLM, the snow depth is limited to Hmax = 1 m liquid water equivalent (LWE), and any additional snow runs off to the ocean. When glacier ice melts, the meltwater remains in place until it refreezes. (In warm parts of the ice sheet, the meltwater does not refreeze but stays in place indefinitely.) In CESM1(CISM), on the other hand, snow exceeding the maximum depth turns to ice, contributing a positive SMB to the ice sheet model. When ice melts, the meltwater is assumed to run off to the ocean, contributing a negative SMB.
In real ice sheets the transition between snow and ice takes place in a firn layer that can be tens of meters thick. Although CLM includes processes important for firn (e.g., refreezing and snow densification), it does not have a true firn model. To better approximate firn, the maximum snow depth Hmax could be increased. We found that for most of the ice sheet, including the ablation zones and the cold interior, the SMB is insensitive to Hmax. However, a greater Hmax results in a more positive SMB in lower-elevation accumulation zones, where a thicker snowpack allows more refreezing and less runoff. Since the CLM snowpack model has not been formally evaluated for large snow depths, we chose the default value Hmax = 1 m, with the caveat that this choice may underestimate refreezing in parts of the ice sheet.
The surface mass balance is usually defined by glaciologists as the total snow accumulation minus the total ablation of ice and snow. In CLM, however, the quantity passed to Glimmer-CISM is the SMB of ice alone. Thus, the snow depth can fluctuate between 0 and 1 m LWE without contributing to the ice sheet SMB. This treatment is adequate for ice sheet evolution on time scales of a few years or longer, but it prevents us from forcing Glimmer-CISM with a realistic SMB on subannual time scales. In future CLM versions, the computed SMB will include changes in both snow and ice mass, in agreement with the usual definition.
c. Coupling of the land and ice sheet models
CESM1(CISM) consists of five physical models (atmosphere, land, ocean, sea ice, and ice sheet) linked by a coupler that merges and regrids fields as needed. The fields sent from the land model to the ice sheet model via the coupler are the surface mass balance (SMB) (computed as described above), surface temperature, and surface elevation in each glacier_mec column of each grid cell. These fields are averaged by CLM over the coupling interval (typically 1 day) before being sent to the coupler. The surface temperature sent to the coupler is the value at the top of the CLM ice column, below any snow. In the case of two-way coupling between CLM and CISM, the surface elevation in each glacier_mec column would change over time. For one-way simulations, however, the surface elevation in each column is fixed at its initial value.
These fields are sent to Glimmer-CISM on the CLM global latitude/longitude grid. They are accumulated and averaged by Glint over the period of a mass balance time step (typically 1 yr) and are downscaled to the ice sheet grid using bilinear interpolation. For each ice sheet grid cell, the SMB and surface temperature are then computed at the local elevation by linear interpolation between the values at the two neighboring elevations. For example, the SMB at 1000 m could be found by averaging the values sent by CLM at 800 and 1200 m. The SMB is applied during the ice sheet thickness calculation, and the surface temperature is used as an upper boundary condition for the vertical temperature calculation.
Communication between the land and ice sheet models is currently one way: Glimmer-CISM receives surface forcing from CLM, but the surface topography and landunit areas in CLM are held fixed. For decade-to-century scale runs, ice sheet geometry changes are modest, and one-way coupling does not lead to large errors. For longer runs with large ice sheet changes, however, two-way coupling will be needed. Much of the infrastructure for two-way coupling is already in place. After each dynamic time step, Glimmer-CISM will send five fields to CLM via the coupler: the ice sheet fractional area and surface elevation, the conductive heat flux at the top surface, and the freshwater fluxes from basal melting and iceberg calving. (The freshwater flux from surface melting is computed directly by CLM.) The fractional area and surface elevation will be used by CLM to update the area of glacier ice in each elevation class; the conductive flux will be applied as a lower boundary condition for the CLM column temperature calculation; and freshwater fluxes will be sent to the ocean model.
3. Model spinup and validation
Here we evaluate the ability of CESM1(CISM) to simulate the Greenland Ice Sheet under preindustrial forcing. We also describe an ensemble exercise to identify parameter sets that combine to produce an optimal GIS, based on an automatic ranking procedure that uses multiple diagnostic measures of model performance.
CESM1(CISM) can be run in several configurations: 1) BG, with all five physical components active; 2) FG, with active atmosphere, land, and ice sheet models; 3) IG, with active land and ice sheet models forced by atmospheric data; and 4) TG, with an active ice sheet model forced by the SMB and surface temperature from a previous BG, FG, or IG simulation. The ice sheet runs described below are TG simulations forced by output from three coupled BG simulations. The BG simulations followed protocols of phase 5 of the Coupled Model Intercomparison Project (CMIP5) (Taylor et al. 2012) and used 1) preindustrial (1850) radiative forcing, 2) twentieth-century (1850–2005) forcing, and 3) twenty-first-century (2005–2100) forcing based on the high-emissions RCP8.5 scenario. As shown in the companion paper by Vizcaíno et al. (2013a), the SMB simulated by CESM for the Greenland Ice Sheet during the late-twentieth century is in good agreement with RACMO2, a state-of-the-art regional climate model (Ettema et al. 2009; van Angelen et al. 2012).
a. Spinup procedure
To simulate future GIS changes, we require a spun-up GIS geometry that reflects the CESM-derived surface mass balance field. It is challenging, however, to spin up a dynamic ice sheet model using SMB forcing from a global climate model. Ideally, ice sheets are driven through at least one glacial cycle before arriving at a preindustrial initial condition, in order to instill an accurate thermal history in the ice (Rogozhina et al. 2011). This is important because internal ice temperatures strongly influence the ice viscosity. But the cost of generating transient paleoclimate forcing fields is proportional to the cost of the driving climate model, and it is infeasible to run CESM over a full glacial cycle.
Use of a surface energy balance model (rather than a positive-degree-day model) to compute surface ablation introduces further difficulties. The classical method of generating climate-model-derived SMB has been to extract annually averaged temperature and precipitation fields from the climate model, add a reconstructed seasonal cycle, and then calculate the SMB using a PDD model forced by these temperature and precipitation fields. With this method one can interpolate between climate “snapshots” (e.g., Last Glacial Maximum and modern) to produce a smooth transient SMB history scaled to paleoproxies such as ice cores. If the same method is applied to future periods, no present-day discontinuity exists. In contrast, snapshot-derived end-member SMB fields constructed using an energy balance model cannot simply be blended because the SMB fields resulting from simple interpolation of cold-end and warm-end members are not physical.
One could in principle use an alternate method (i.e., a simple ice-core-derived climate parameterization or the PDD-derived snapshot approach) to spin up the ice sheet model through a glacial cycle, then switch to energy balance forcing for simulations of the future climate (e.g., Vizcaíno et al. 2010). However, this method yields a discontinuity in forcing at the time of switchover. The discontinuity results in perturbations to ice geometry that contribute to a nonphysical transient component of ice sheet evolution. Also, this method does not allow for explicit validation of the preindustrial CESM-derived SMB since the spun-up ice sheet does not depend on CESM forcing.
For the simulations described below, we chose to spin up Glimmer-CISM in TG mode with SMB and surface temperature from a coupled CESM simulation with preindustrial (1850) forcing. This BG1850 simulation used active atmosphere, land, ocean, and sea ice components at ~1° resolution. The ice sheet model was run at 5-km resolution, with the ice sheet SMB computed as described in section 2b. The GIS was initialized to modern geometry using SeaRISE datasets (http://websrv.cs.umt.edu/isis/index.php/SeaRISE_Assessment) and was not allowed to evolve. For the last 75 years of the simulation, we stored the annual-mean SMB and surface temperature in each elevation class of each glaciated grid cell. We then applied this forcing to CISM in a repeat cycle over many model millennia, during which the ice sheet geometry and temperature evolved toward steady state under CESM-simulated interannual and interdecadal variability.
This method neglects the thermal signal of previous glacial periods, which may affect the model optimization exercise (section 3b). It does, however, allow for generation of smooth, continuous forcing between historical and future periods. In addition, we suggest below that internal temperatures likely have little effect on twenty-first-century mass loss predictions generated with this model using optimized geometry because of the lack of dramatic differences in the near-term ice dynamic response for different ice sheet model configurations. Thus, while a transient spinup procedure is a logical next step in model development and may affect the choice of parameters derived in the optimization exercise, we suggest that a different procedure would not strongly affect the results in section 4.
b. Ensemble-based model optimization
To generate a preindustrial spun-up GIS with optimized geometry, we used an ensemble approach rather than manual tuning. Others have recently used this methodology to optimize ice sheet simulations and explore model parametric uncertainty (Stone et al. 2010; Applegate et al. 2012), and we follow their general approach here. First, we chose three parameters that have been suggested to provide important controls on ice sheet model behavior: the flow enhancement factor (f), the coefficient of the basal traction equation (B), and the geothermal heat flux (G). The scalar f controls the ice viscosity (with large f implying lower viscosity) and parameterizes the effects of impurities and anisotropy on ice flow. The coefficient B relates the basal sliding velocity to the driving stress and thus parameterizes subgrid-scale processes such as bed roughness and basal water pressure. We assumed a linear relation between sliding velocity and driving stress where the bed is thawed (with large B implying faster sliding for a given stress) and no sliding over a frozen bed. The flux G influences where the bed is frozen or thawed. As in previous studies, we simplified the analysis by using a single value of each parameter for the entire ice sheet, even though these parameters are known to be spatially heterogeneous. The parameter ranges (Table 1) were chosen based on preliminary CISM simulations and other studies (e.g., Stone et al. 2010; Applegate et al. 2012). Within the range of each parameter, we assumed that all possible values were equally likely: in other words, we did not ascribe an a priori probability distribution function to parameter sets. This is because our main goal was to identify optimal model configurations by sampling the entire range, rather than to explore parametric uncertainty.
Based on these ranges, we generated 100 parameter sets using Latin hypercube sampling (LHS) and used each parameter set in a 21 000-yr Glimmer-CISM simulation. The model was run on a 5-km Greenland grid with a dynamic time step of 0.1 yr, decreased to 0.05 yr for the last 9000 years of the simulation to damp numerical waves in the Jakobshavn region. Surface forcing for all simulations was identical and was derived from the BG1850 coupled run. Thus the SMB was not a free parameter in the optimization process. Rather, we aimed to find the ice sheet parameters resulting in the most realistic GIS simulation for the given SMB forcing. The number and length of runs were determined by the computational expense required for the ice sheet model to reach quasi equilibrium while still adequately sampling parameter space. Internal temperatures were initialized with a linear vertical profile; we prescribed observed modern temperatures at the surface and a value 2°C below pressure melting point at the bed. While internal temperatures were still drifting residually at the end of the simulation, the 21 000-yr simulation resulted in geometric near equilibrium (Fig. 2a) for almost all ensemble members, indicating that temperature drift had ceased to noticeably affect geometric evolution. The final state of each simulation was compared to an observational dataset consisting of the Bamber et al. (2001) Greenland dataset, augmented with high-resolution ice thickness data around Jakobshavn Isbrae (Gogineni 2012). The approximate areas of glaciers and ice caps (GIC) peripheral to the Greenland Ice Sheet were calculated from the Randolph Glacier Inventory (Arendt et al. 2012), with GIC volume estimated using volume–area scaling relationships (Bahr et al. 1997). The modeled GIS was evaluated against these data based on five diagnostics: ice volume, ice area, rms error (RMSE) of ice thickness, maximum (summit) elevation, and horizontal summit offset. Total observed area and volume were derived from the Bamber and Randolph Glacier Inventory data, while observed surface elevations used in RMSE and summit diagnostics were derived solely from the Bamber et al. (2001) dataset. Although absolute RMSE values for ice thickness may contain artifacts related to autocorrelation, we believe that relative RMSE errors are still a reasonable measure of model performance. To determine the ensemble member (and thus parameter combination) that best fit these diagnostics, we developed an objective ranking algorithm that automatically sorted ensemble members according to their worst-ranked diagnostic. This “equally weighted worst diagnostic” approach punishes simulations that do very well for one diagnostic at the expense of others, while rewarding simulations that do reasonably well for all diagnostics.
Given the simplified ice sheet physics (e.g., the shallow-ice approximation, the linear sliding law, and the spatially uniform model parameters), we would expect significant errors in equilibrium geometry even with a perfect SMB. We show below, however, that the SMB is an important first-order control. A poor SMB would yield an unrealistic ice sheet geometry even with sophisticated physics. CESM gives a good SMB simulation for most of the GIS (Vizcaíno et al. 2013a), but the biases that exist in originally ice-free regions, especially near the northern and eastern margins, result in an ice sheet that is too thick and extensive for all parameter choices. Thus, in addition to highlighting ensemble members that best reproduce the observed ice sheet, the optimization exercise identifies CESM biases that are independent of the ice sheet model.
c. Common features of spinup simulations
The 100 ensemble members showed some consistent responses to CESM-derived forcing. Most notably, the total ice area and volume were always overestimated (Fig. 2b). The degree of overestimation ranged from 11%–33% for volume and was roughly 18% for ice area; no ensemble member underestimated the volume or area. The excess ice area and volume resulted primarily from anomalous ice growth around the margins (Fig. 3). This suggests that 1) dynamic ice expansion was excessive and/or 2) the SMB in regions of margin advance was too positive, driving in situ ice growth where no present-day ice is observed. To evaluate the importance of 2), we examined the SMB for an arbitrary ensemble member after one year of model integration, when the ice sheet was close to modern geometry and marginal ice growth had yet to occur. The most extensive regions that are ice free at present but were consistently ice covered in the spun-up model were Peary Land (section aa′ in Fig. 3), the southeast (Fig. 3, east side of section cc′), and the Geikie Plateau margins. In these regions the simulated SMB at the start of the simulation was slightly positive (Fig. 4), resulting in perennial snow cover and subsequent in situ ice growth in these initially ice-free locations. As the spinups progressed, the new marginal ice merged with the main ice sheet.
In some modeled regions with excessive ice extent, there is indeed extensive present-day glacier or ice cap cover that is not captured in the Bamber et al. (2001) dataset, but is documented by Arendt et al. (2012). It is therefore not unrealistic to see sporadically positive SMB values in these regions. However, the presence of a broadly uniform positive SMB combined with the general absence of net ablation zones even at sea level indicates that SMB bias in these regions (even if small) largely determines the final GIS geometry. Attributing and correcting these biases in the CESM Arctic climate simulation will be a priority for future work.
In contrast to the marginal regions that experienced anomalous in situ ice growth, there is a smaller region of excess ice growth on the southwest margin (section cc′ in Fig. 3, west side). This growth can be attributed primarily to dynamic ice advance, given a zero SMB over bare land adjacent to an established ablation zone at the start of the spinup. (Zero SMB over bare land indicates that seasonal snow does not survive the summer.) This dynamic advance likely resulted from anomalously weak ablation or excessive inland accumulation in the region. In both the in situ and dynamic excess ice cases, buttressing of interior ice and instigation of ice height–SMB feedbacks reinforced marginal ice growth as simulations proceeded.
Although small positive SMB values occurred over parts of the GIS margin that are currently ice free, the SMB on the whole over the existing ice sheet is very realistic, as discussed by Vizcaíno et al. (2013a). In regions where excess marginal ice did not develop, simulated ice sheet geometry compared well with observed geometry, with differences among ensemble members due to varying ice sheet parameters (e.g., cross sections aa′ and bb′ in Fig 3). In addition, despite the simple shallow-ice approximation (SIA) dynamics, the 5-km ice sheet grid generally allowed focused ice flow into outlet glaciers, in good agreement with the spatial pattern of balance velocities based on observations (Fig. 5). This was because the basal topographic dataset was able to resolve outlet glacier valleys, allowing thicker ice and faster flow to develop locally. The simulated flow in these regions, however, was dominated by vertical shear, with little basal sliding. As a result, the maximum ice speeds in the modeled outlet glaciers were consistently lower than observed. A higher-order flow model and spatially varying basal traction field (e.g., Price et al. 2011) would likely be needed to simulate ice speeds accurately in outlet glaciers.
d. Model ranking
After spinup, each ensemble member was ranked according to the diagnostics described in section 3b. Figure 6 illustrates the rankings for each diagnostic. As noted above, area and volume were consistently overestimated, especially along the northern and eastern margins. The RMSE of surface height ranged from 240 to 380 m with significant errors in the regions of excess marginal ice growth. Summit elevations were, on average, slightly lower than observed. That is, most of the excessive ice occurred at the margins and not in the interior. Finally, the horizontal offset of the summit elevation displayed an uneven bimodal structure: most simulations had an offset of 75–100 km, while a few poorly ranked models had horizontal offsets of nearly 500 km. The large offsets in these poorly ranked simulations resulted from lowering of the main GIS dome to elevations lower than in the southeast mountains, highlighting these simulations as having an unrealistic geometry.
To explore the relative importance of each Latin-hypercube-sampled parameter on individual diagnostics, the parameter and diagnostic values were plotted against each other (Fig. 7). The ice flow enhancement factor f played a large role in model performance across diagnostics. The basal traction parameter B was generally of secondary importance. The dependence of the diagnostics on B was greatest at higher values of B (i.e., lower traction and faster sliding); at lower values, small changes in B were overwhelmed by variations in f. Furthermore, the effects of changes in B (at low B) were concentrated at the margins, with little impact on large-scale ice sheet geometry. At least in the range tested here, G had little influence on any diagnostics, in agreement with earlier studies (Stone et al. 2010; Applegate et al. 2012). A spatially varying G field, however, could affect the GIS geometry in ways that are not directly testable here (e.g., through its influence on subglacial hydrology and basal sliding).
We used the equally weighted worst-diagnostic ranking approach (section 3b) to identify simulations that performed well in all diagnostic categories. The top five simulations are shown in Fig. 7, and their simulated GIS geometries are shown in Fig. 8. The best simulations had a flow enhancement factor f clustered between 2.9 and 3.2. With one exception, the optimal models had B values of less than ~10−5 m yr−1 Pa−1. The nonlinear horizontal axis in Fig. 7 reflects our use of the exponent of the basal traction parameter (rather than the parameter itself) in the Latin hypercube sampling. The axis thus understates this clustering. If a linear axis were used, the majority of the optimal values would collapse to a tight cluster at low B. Simulations with high B tended to have large vertical and horizontal summit offsets, reflecting a “slumped” profile from excessive sliding. This suggests that a spatially varying B field, perhaps obtained through inverse modeling, would better capture ice sheet geometry through lower B in the interior and higher B in outlet glaciers. The top five simulations have a wide range of G, showing that the choice of G within the prescribed bounds has little effect on the ranking.
Internal temperatures (which are likely too warm because of the spinup procedure) may have affected the choice of best parameter sets. In particular, a colder ice sheet would require a larger value of f to avoid an overly thick and viscous ice sheet. This agrees with the observation that our best values for f are relatively low compared to empirical estimates (Cuffey and Paterson 2010).
In no cases did “optimal” models score first in any one diagnostic category. We consider this a success for the following reasons. Figure 9 compares the overall top-ranking simulation to the simulation that ranks best for ice volume alone. In the latter, marginal ice growth due to SMB anomalies is compensated by large values of f and B, resulting in fast outflow from the interior, unrealistic geometry, and low rankings for horizontal and vertical summit offset. The use of multiple diagnostics identifies this as a poor simulation and minimizes the effect of SMB anomalies on the choice of optimal parameter sets. However, it does not entirely remove the effect of SMB anomalies; each of the best five simulations underestimated the summit elevation to compensate for excess marginal growth by increasing either f or B.
The optimal f, B, and G parameter sets were potentially affected by biases from the CESM-generated surface mass balance, the spinup procedure, and the use of nonspatially varying (scalar) parameters. If these biases were reduced, the optimization exercise would likely yield a different set of best parameter values.
4. Future GIS evolution under CESM forcing
An important application of CESM1(CISM) is to estimate sea level rise resulting from GIS evolution over the historical period (1850–2005) and through the twenty-first century. To provide an initial estimate and explore the effect of ice sheet parameter choices on sea level projections, we forced Glimmer-CISM with the transient SMB and surface temperature from two coupled (BG) CESM simulations: a historical twentieth-century simulation (1850–2005) followed by a twenty-first-century representative concentration pathway (RCP8.5) simulation (2005–2100). We then evaluated changes in GIS geometry resulting from SMB changes and the relative role of ice dynamics and SMB on the GIS sea level contribution.
a. Historical volume change
Starting from the ice sheet state at the end of the preindustrial spinup, we ran each of the five top-ranking ensemble members forward from 1850–2100. The climatological SMB for the simulated GIS during the early historical (1850–80), modern (1970–2000), and future (2070–2100) periods, averaged over these ensemble members, is shown in Fig. 10. Time series of area-integrated SMB, ice discharge, and volume evolution are shown in Figs. 11a–c. Volume change is given in units of sea level equivalent, where 1 mm of sea level corresponds to about 360 Gt of water. There is little variation in trend among the five top ensemble members. During 1850–2005, a gradual decrease in ice volume accounts for 16 mm of sea level rise or 7% of the 1880–2009 global average SLR reconstructed by Church and White (2010). The spatial pattern of ice thickness change resulting from the historical SMB (see the 1970–2000 average SMB on Fig. 11a) is shown in Fig. 12a. Ice loss occurred mainly through ablation near the margins; this loss was slightly compensated by elevation increases in the interior due to increased accumulation (Vizcaíno et al. 2013a).
b. Future volume change
During 2005–50 the rate of sea level rise from GIS mass loss was 1.9 mm decade−1 as a result of a more negative surface mass balance, as shown in Fig. 11a. The SMB declined more sharply after 2050, accompanied by an increase in SMB variability. After 2070 the 20-yr moving average SMB became negative, implying long-term decay of the GIS even if discharge were reduced to zero. By the last three decades of the simulation, the GIS contribution to SLR increased to 13.3 mm decade−1. The average total GIS contribution to SLR for the top-ranking ensemble members was ~76 mm between 1850 and 2100, with 60 mm occurring after 2005. We emphasize that this is a single climate model instance; other instances would have different interdecadal variability superposed on the long-term trends.
The spatial pattern of twenty-first century ice loss was qualitatively similar to that of 1850–2005 but greater in magnitude (Fig. 12b). The ice thinned around the margins, in some regions by more than 400 m (mainly where the ice sheet retreated, leaving bare land). In the interior, thickening owing to increased accumulation raised the elevation by 5–10 m over large areas of the ice sheet. Margin retreat without summit lowering implies increasing surface gradients along the flanks. During a longer simulation, steeper surface slopes would drive faster ice flow from the interior, causing dynamic lowering that would likely overwhelm the effect of increased accumulation.
Figure 11b shows a decreasing trend in ice discharge to the ocean during 1850–2100. This decrease could result from one or more of three mechanisms: 1) retreat of the ice margin to a grounded position, reducing the area of marine-terminating ice; 2) ice thinning or reduction of the surface gradient, resulting in slower flow near marine margins (due to the nonlinear relationship between ice speed and thickness/surface slope in the shallow-ice approximation); and 3) increased ablation, resulting in less ice available for discharge. Mechanisms 2 and 3 are closely related because increased ablation will reduce discharge across the calving front via direct thinning as well as decreased velocity because of mechanism 2. To explore the role of these mechanisms in decreasing discharge, we analyzed the ice geometry and fluxes near the coastal margins. At the boundary of the spun-up 1850 ice sheet, 58% of grid cells bordered the ocean. From 1850–2100 this total decreased by a modest 6%, showing that the transition of the GIS to a dominantly land-terminating ice sheet would require several more centuries of simulation. However, ice flux into these ocean-terminating grid cells decreased by about 30%, and the upstream ablation zone expanded significantly, showing that both mechanisms 2 and 3 were major factors in lowering the discharge.
Figure 11c illustrates the sea level contributions from SMB changes and dynamic feedbacks. The middle black/orange line shows the cumulative SLR anomaly from GIS mass loss, bracketed by a blue line showing the cumulative SLR anomalies due to SMB alone and a red line showing the difference between the black/orange and blue lines. The residual signal (the red line) shows the cumulative effect of changes in dynamic discharge. The decrease in discharge attributed to ice thinning is a large negative feedback, offsetting ~50% of the mass loss due to a more negative SMB.
The magnitude of this negative feedback should be interpreted with caution. The anomalous advance of the ice sheet margin during spinup gives an excessive fraction of marine-terminating ice sitting low in the future ablation zone, likely influencing the size of the feedback. Moreover, the transient response of outlet glaciers could differ with a higher-order ice flow model or with a more advanced sliding law that allows positive feedbacks on ice flow. Our results are consistent, however, with a recent study using a higher-order model (Gillet-Chaulet et al. 2012), which found a large decrease in discharge in response to twenty-first-century warming and melting (leading in turn to thinning and reduced driving stresses).
Our model does not capture some potential positive dynamic feedbacks. Since there is no direct coupling between the ice sheet and the ocean, acceleration and thinning of outlet glaciers owing to ocean warming (e.g., Holland et al. 2008; van den Broeke et al. 2009) were not simulated. Also, potential ice flow changes attributed to evolving subglacial hydrology (e.g., Hoffman et al. 2011; Sundal et al. 2011) were not resolved.
c. Role of ice sheet parameters in SLR projections
The SLR projections in section 4b were based on the five top-ranking simulations. To explore the dependence of sea level rise projections on ice sheet model parameters, we repeated the analysis using the five worst-ranking ensemble members, whose spun-up geometry was unrealistic as measured by one or more diagnostics. When integrated from 1850 to 2100, these ensemble members gave net SLR ranging from 60 to 87 mm, compared to 75–78 mm for the tightly clustered top-ranking members. Notably, the surface mass balance among the five low-ranking ensemble members diverged after 2000, but there was much less growth in the related spread of discharge rates. The models with the largest negative changes in SMB were those with relatively low surface elevation because of low viscosity combined with low basal traction. These models had anomalously large increases in ablation area for a given lifting of the equilibrium line altitude (ELA). This result suggests that the SLR response is sensitive to the initial ice sheet geometry but is less sensitive to the choice of dynamic parameters, provided that these parameters yield a realistic spun-up geometry. Thus, choosing dynamic parameters that reproduced the modern GIS geometry as closely as possible was important for near-term SLR projections, not because the parameter values themselves yielded major dynamic changes but rather because more accurate geometry resulted in more realistic SMB sensitivity to a higher ELA. On longer time scales, ice sheet parameter choices would have a growing effect on the dynamic response.
5. Conclusions and future development
We implemented a dynamic ice sheet model, Glimmer-CISM, in a coupled global climate model, CESM. The surface mass balance (SMB) is computed in multiple elevation classes per grid cell in the land model and then downscaled to the ice sheet model. When Glimmer-CISM is forced with the SMB from a coupled preindustrial CESM run, the simulated steady-state Greenland Ice Sheet is in broad agreement with the observed geometry, although it is too thick and extensive along the northern and eastern margins. The spatial pattern of simulated velocities is similar to observationally based balance velocities over most of the ice sheet, despite the use of the shallow-ice approximation, but maximum ice speeds in outlet glaciers are too low.
To explore the dependence of the steady-state geometry on ice sheet parameters, we carried out a 100-member ensemble of 21 000-yr CISM simulations. Each simulation was forced with the same preindustrial SMB from CESM but had different values of ice sheet parameters associated with internal viscosity, basal sliding, and geothermal heat flux. We used an objective, multidiagnostic ranking scheme to identify simulations that gave the best fit to observed GIS geometry. The ice flow enhancement factor f, which controls internal viscosity, played a leading role in determining the steady-state GIS geometry; values of f for the top-ranking simulations were closely clustered from 2.9 to 3.2. Low values of basal sliding also tended to produce more realistic ice sheet geometry. Geothermal heat flux was of little importance over the range tested. All simulations yielded a steady-state ice sheet with greater area and volume than the observed GIS, mainly because of an anomalous positive SMB in some regions that are currently ice free. Optimal parameter settings would likely change if these SMB biases were reduced.
Following spinup, we forced the five top-ranking ensemble members with the SMB from coupled CESM simulations of the twentieth century (1850–2005) and twenty-first century (2005–2100). The simulated sea level rise associated with GIS mass loss was 16 mm during 1850–2005 (with historical radiative forcing) and 60 mm during 2005–2100 (with RCP8.5 forcing). The RCP8.5 projection is within the range of other recent sea level assessments using different methodologies (e.g., Bindschadler et al. 2013; Shannon et al. 2013). This mass loss was driven primarily by increased ablation near the margins. Sea level rise from a decreasing SMB was offset by reduced ice discharge to the ocean, caused mainly by thinning and slowing of ice near coastal margins. Sea level rise (SLR) projections are sensitive to the initial GIS geometry, and therefore it is important to choose model parameters that give an accurate spun-up geometry.
The coupled CESM1(CISM) model requires further development to reduce biases and increase realism. Desired model improvements can be divided into four categories: 1) a more sophisticated dynamic ice sheet model, 2) more accurate downscaling of atmospheric fields to the land surface, 3) more realistic treatment of snow and ice physics in CLM, and 4) two-way coupling between CISM and other CESM components.
The version of Glimmer-CISM in CESM uses the shallow-ice approximation. Developmental versions of CISM include several higher-order velocity solvers (Price et al. 2011; Perego et al. 2012; Cornford et al. 2013) that are valid for fast flow and will be included in future versions of CESM. Meanwhile, we are testing nonlinear sliding laws (e.g., Schoof 2005) that incorporate plastic flow and interactions with basal hydrology, in order to obtain better agreement with observations of fast-flowing ice. We will also implement improved schemes of iceberg calving (e.g., Levermann et al. 2012). The ice sheet model will then be able to support simulations not only of the GIS but also of the Antarctic Ice Sheet, much of which is marine based and bordered by large ice shelves.
The current downscaling scheme converts the temperature and specific humidity at the mean gridcell elevation to values for each elevation class, assuming a fixed globally uniform lapse rate. Surface radiation and precipitation fields are not downscaled to elevation classes. We are now testing schemes to downscale the downward longwave radiation based on the temperature lapse rate and to convert precipitation between rain and snow based on the downscaled surface temperature. (These schemes require some additional code changes to conserve energy.) Spatially and temporally varying lapse rates (e.g., Jarosch et al. 2012) also could improve the accuracy of the SMB. In the longer term we will run CESM using a version of the Community Atmosphere Model (CAM) with regional grid refinement over ice sheets so as to improve the simulation of orographically forced precipitation and to provide more accurate atmospheric fields prior to downscaling.
The Community Land Model snowpack model includes detailed treatments of snow metamorphosis, meltwater percolation, and refreezing. Bare ice, however, is treated crudely, with a prescribed albedo for visible and near-IR radiation bands. Observed ice albedos are sensitive to the presence of surface meltwater and impurities, and thus a more realistic treatment of radiative transfer in ice could improve the simulation of shortwave absorption and ice ablation. Also, the snowpack model was developed primarily for seasonal snow with a depth of no more than a few meters. Densification and refreezing could be simulated more accurately with a firn model that resolves a snow/ice transition depth of tens of meters.
Finally, the current model sends the ice-sheet surface mass balance from CLM to Glimmer-CISM, but the evolving ice sheet geometry is not returned to CLM. To support simulations on century-to-millenial time scales, we are now modifying CLM to support two-way coupling with ice sheets. CLM surface topography and land surface types will evolve, with glaciated landunits converted to vegetation and vice versa as ice sheets retreat and advance. Work is ongoing to allow CAM atmospheric flow to respond dynamically to changes in ice sheet topography. We also plan to couple Glimmer-CISM to the CESM ocean model, with the ultimate goal of simulating interactions of ice sheets with all other components of the climate system.
We thank Gail Gutowski and Charles Jackson for valuable interactions and advice during model testing. We also thank Janneke Ettema, Jonathan Gregory, Magnus Hagdorn, Matthew Hoffman, David Holland, Jesse Johnson, Tony Payne, Stephen Price, Jeff Ridley, Ian Rutt, and Michiel van den Broeke for useful discussions. Support for WHL was provided by the Scientific Discovery through Advanced Computing (SciDAC) project funded by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Biological and Environmental Research (BER). JGF was supported by the Regional Arctic System Modeling project funded by BER and the National Science Foundation's Office of Polar Programs. The Los Alamos National Laboratory is operated by the DOE National Nuclear Security Administration under Contract DE-AC52-06NA25396. MV was supported by a Marie Curie Incoming International Fellowship at Utrecht University. Support was also provided by the National Science Foundation through Awards ATM-0917755 (for MV), ANT-1103686 (for WJS), and a Small Grant for Exploratory Research (for WHL). The CESM project is supported by the National Science Foundation and the Office of Science (BER) of the U.S. Department of Energy.
Computing resources were provided by the Climate Simulation Laboratory at the NCAR Computational and Information Systems Laboratory (CISL), sponsored by the National Science Foundation and other agencies. The CMIP5 simulations used in this research were enabled by CISL compute and storage resources. Bluefire, a 4064-processor IBM Power6 resource with a peak of 77 TeraFLOPS, provided more than 7.5 million computing hours, the GLADE high-speed disk resources provided 0.4 PetaBytes of dedicated disk, and CISL's 12-PB HPSS archive provided over 1 PetaByte of storage in support of this research project. In addition, this work used resources of the Oak Ridge Leadership Computing Facility, located in the National Center for Computational Sciences at Oak Ridge National Laboratory, which is supported by the Office of Science of the Department of Energy under Contract DE-AC05-00OR22725. We acknowledge the use of data and/or data products from CReSIS generated with support from NSF Grant ANT-0424589 and NASA Grant NNX10AT68G.
The National Center for Atmospheric Research is sponsored by the National Science Foundation.
This article is included in the CESM1 Special Collection.