Significance This study integrates an economic model of climate change with a small structural model of the Greenland ice sheet (GIS). As such, it provides a methodology for incorporating large earth system changes into standard economic cost–benefit or damage-limiting analyses. It finds that adding the GIS has only a small effect on the social cost of carbon (SCC) because melting is slow and damages are far in the future.

Abstract Concerns about the impact on large-scale earth systems have taken center stage in the scientific and economic analysis of climate change. The present study analyzes the economic impact of a potential disintegration of the Greenland ice sheet (GIS). The study introduces an approach that combines long-run economic growth models, climate models, and reduced-form GIS models. The study demonstrates that social cost–benefit analysis and damage-limiting strategies can be usefully extended to illuminate issues with major long-term consequences, as well as concerns such as potential tipping points, irreversibility, and hysteresis. A key finding is that, under a wide range of assumptions, the risk of GIS disintegration makes a small contribution to the optimal stringency of current policy or to the overall social cost of climate change. It finds that the cost of GIS disintegration adds less than 5% to the social cost of carbon (SCC) under alternative discount rates and estimates of the GIS dynamics.

The future of the mammoth Greenland ice sheet (GIS) is one of the largest and most complicated issues facing environmental policy in the coming years. Complete disintegration of the GIS would raise the level of the oceans by more than 7 m, inundating many of the world’s major human settlements. Paleoclimatic findings, as well as ice sheet modeling, indicate that the current trajectory of global temperatures would lead to nearly complete disappearance of the GIS over the coming millennia. The critical questions are, how fast will the ice sheet decline, and what can be done to stop the disintegration and resulting inundation?

The present study examines economic aspects of the disintegration of the GIS by incorporating a small reduced-form model of the GIS into the Dynamic Integrated model of Climate and the Economy (DICE) of the economics of climate change. Studies find that a rise in temperature a few degrees above the current levels will lead to a nearly total ice sheet loss. The Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (1) and the IPCC 1.5 °C report (2) cite several studies that suggest nearly complete melting will occur at thresholds as low as 1 °C and as high as 4 °C of warming. One prominent study, by Robinson et al. (3), estimates “that the warming threshold leading to a monostable, essentially ice-free state is in the range of 0.8–3.2 °C, with a best estimate of 1.6 °C.” (Note that studies use different reference years in defining temperature increases. For present purposes, I define the “base temperature” as that global temperature associated with a stable GIS. The GIS was stable in the 1961–1990 period, so T = 0 is associated with global mean temperature in that period. The standard definition of “pre-historic” temperature adopted by ref. 2 is for the 1850–1900 period. By that standard, the GIS-stable “base temperature” adopted here is about 0.3 °C above prehistoric temperatures.)

Further complicating the question is uncertainty about whether there are multiple sets of equilibria of temperature and GIS volume, but the best modeling evidence suggests multiple equilibria with hysteresis (see refs. 3 and 4 as discussed below).

An additional unresolved factor is the dynamics of disintegration and rebuilding between equilibria. Modeling studies indicate that the path of melting is slow, with the central estimate being that, at a 6 °C global warming, the GIS would lose 10% of its volume in four to five centuries. The exact dynamics vary widely among alternative models.

The current study develops a model of GIS equilibrium and dynamics that is based on current studies but sufficiently compact to integrate fully into an economic model. The result is the DICE-GIS model, which includes the standard components of the DICE-2016R2 integrated assessment model. Based on this combined model, the study then examines baseline (no-climate-policy) and optimal climate policies along with different constraints, parameters, and discount rates.

Here are the major results of the study. First, the study shows that integrated economic–geophysical modeling is a promising approach to policy analysis of major earth system changes. While the study applies only to the GIS, the same method can in principle be used for other major and potentially catastrophic changes, such as those relating to the Antarctic ice sheet or changes in the North Atlantic thermohaline circulation. Integrated modeling will help illuminate alternative policies and results as well as the design of global strategies to prevent catastrophic changes.

Second, the study finds that a baseline or no-policy path will lead to the gradual melting of the GIS over the coming millennium. The GIS reaches the Robinson upper tipping point of 80% (explained below) in 500 y. Once past the tipping point, disintegration is difficult to reverse. On the other hand, strong climate policy in the optimal run can stop the GIS decline well short of complete disintegration or critical tipping points. Put differently, under a wide range of assumptions, policies based on the social cost of carbon (SCC) used by governments today are likely to slow warming sufficiently to reduce the risk of irreversible GIS disintegration.

Third, there are two alternative approaches to including the GIS in policy studies: either through including a damage function or through putting a volumetric limit on the decline of the GIS. These give roughly the same policies for the GIS dynamics for most discount rates and central estimated melt rates. The study also shows the usefulness of cost–benefit analyses as a method for developing estimates of the SCC, the optimal carbon price, and optimal/efficient emissions pathways.

Fourth, the most important single indicator of the strength of current climate policy is the SCC. Therefore, a useful way of understanding the impact of GIS disintegration on climate policy is to estimate how much adding GIS damages or GIS volumetric constraints changes the SCC. This study demonstrates that, under a very wide range of assumptions, the risk of GIS disintegration—although a major change in the earth system—would make a small further contribution to the overall SCC or to the overall cost of climate change. The increment to the SCC is near zero at moderate discount rates and as high as 5% of the total SCC at very low discount rates and high melt rates. At the discount rate used by the US government, the addition of GIS damages to the SCC is essentially zero. The intuition behind this result is that the timescale of damages associated with GIS melting is much longer than those associated with non-GIS damages. For example, impacts on agriculture are determined largely by contemporaneous climate changes, while GIS melting is determined by climate only with a long lag.

Fifth, the study considers four alternative approaches to dynamics: linear, nonlinear monotonic, irreversible, and hysteretic. The basic finding is that to a first approximation the results are identical for all approaches. In other words, for the GIS, the optimal policies are essentially independent of whether the dynamics are simple (linear) or complex (including hysteretic with multiple equilibria).

Sixth, the present study includes only one of several potential large-scale earth system changes that threaten key physical, biological, chemical, and ecological systems on which human and natural systems depend. Other important systems include the Antarctic ice sheet, Atlantic meridional overturning circulation changes, ocean carbonization, monsoons, tropical cyclones, and forests. The present study may overestimate or underestimate the economic impacts of the GIS depending upon the nonlinear interactions between the GIS and these other earth system changes. For example, if the economic impact of sea level rise is highly convex with respect to the rate or level of the rise, then the impacts here will be underestimates because the calculated marginal impacts are smaller than the true ones that include other impacts of warming. This shortcoming is an important qualification but must await studies of other systems to be evaluated.

Finally, the consideration of geoengineering options leads to an important implication. Simulations with DICE-GIS and large ice sheet models (ISMs) indicate that there is a sharp asymmetry in the response of ice sheet changes to high and low temperatures. A policy that reduces global and GIS temperature to preindustrial levels produces a very slow rebuilding of the GIS compared with the pace of melting at high temperatures. From an economic and policy perspective, the implication is that GIS disintegration should be viewed as an irreversible process. While the GIS may eventually rebuild to its current volume if temperature declines to preindustrial levels, the rate of rebuilding is so slow that the damage cannot be undone within the time perspective of climate policy and human settlements.

Modeling Complex Dynamic Systems Concerns about the impact of climate change on large-scale and unmanageable earth systems have taken center stage in the scientific and economic analysis of climate change. Continued warming threatens to push large-scale earth systems beyond tipping points (2, 5, 6). In the present context, a tipping point is an unstable equilibrium or one from which small shocks will lead to significantly different long-run equilibria. This is also sometimes expressed as saying that small shocks can have large consequences for a system. Important examples for earth systems are the potential collapse of the North Atlantic thermohaline circulation, the dieback of the Amazon rainforest, species losses, and melting of the giant ice sheets. In modeling such complex dynamic systems, it will be important to examine different kinds of equilibria, beginning with the basic structure in the next section. A key question for the present study is the structure of the equilibria of the GIS (or analogous dynamic systems). Three important types are systems that are reversible, hysteretic, and irreversible. Here are simple definitions of these processes: A reversible system is one with no memory, as with a stick that bends and then returns to its original position. An irreversible or highly asymmetrical system is one that breaks or changes to a new state once a threshold is passed, as with a stick that is broken when bent too far. A hysteretic system is one with path dependence or memory of its history. Here, an example would be the consequence of an abrupt climatic event. With a given climate, certain species (such as dinosaurs) would thrive, whereas, after a sharp climatic change (such as a sharp cooling for a few centuries), an entirely new ecosystem might evolve when the climate returned to its original state. Focusing on ice sheets, begin with an equilibrium relationship between temperature (T*) and ice sheet volume (V*): T * = f ( V * ) . [1]The function is written in this form because there is a unique temperature that is associated with every equilibrium GIS volume. However, there may be multiple GIS volumes associated with a given temperature. In the initial specification used in this study, the inverse function is one-to-one. However, in specifications investigated later, such as ones displaying hysteresis, the inverse function does not hold uniquely (that is, there may be multiple equilibrium volumes associated with a single temperature). Fig. 1 shows three alternative versions of Eq. 1. Fig. 1A shows a completely reversible dynamic system, where the ice sheet marches down the f(V*) curve in a warming world, and then marches back up along the same curve as temperatures fall. Fig. 1B is the hysteresis diagram, which will be examined in a subsequent section. The irreversible case in Fig. 1C can be interpreted as an extreme example of the hysteresis curve, where the lower branch has its lower inflection point B at an extremely low temperature. Similarly, the reversible case is at the other extreme of hysteresis where the two branches in Fig. 1B collapse into one branch. Fig. 1. Alternative specifications of GIS equilibrium. A is reversible; B displays hysteresis; C is effectively irreversible because rebuilding requires ice age conditions. The GIS appears to fall into the hysteretic category in Fig. 1B. The hysteretic case presents difficulties of interpretation and modeling. The section between points A and B will not be observed in equilibrium as they are dynamically unstable. Moreover, for temperatures between A and B, there are multiple stable equilibria. The (T*, V*) equilibria will move down the upper branch in Fig. 1B as the ice sheet melts, but the equilibria will follow the lower branch up and to the left as temperature rises and the GIS rebuilds. The plan of the paper is the following. It begins with a discussion of the mathematical and empirical structure of the GIS. This is followed by modeling details and results. The subsequent sections analyze the impact of alternative equilibrium structures. The final section presents reservations and qualifications.

Further Analysis of the GIS Existing Studies. The scientific literature on the impacts of major or catastrophic changes in earth systems is vast. In the scientific domain, the 2013 IPCC (1) reviewed several potential major “abrupt” changes but found “low confidence and little consensus on the likelihood of such events over the 21st century.” The 2018 IPCC report (2) reviewed findings since ref. 1 and found a temperature increase of 1.5–2 °C may be regarded as representing moderate risk for large ice sheets (these findings are reviewed in detail below). There is also a growing set of studies integrating economics and large-scale earth systems, such as collapse of the Atlantic circulation and ice sheet melting (7⇓⇓–10). One example is the development of a calibrated model (“SIMPLE”) similar to the one used here to test the impacts of geoengineering (10). However, integrated economic–climate models of tipping points and catastrophes have been schematic and have generally not relied on realistic geophysical models. This section begins with a review of current physical GIS models and what they suggest about the dynamic structure of giant ice sheets. It then develops a manageable dynamic model of the GIS and explains how to include that in the DICE integrated assessment model. Physical Models of the GIS. A brief description of the GIS may be useful for nonspecialists. Greenland is the world’s largest island, with an area of 2.17 million km2 or about five times the size of California. The ice sheet covers 1.7 million km2, or about 80% of the area, with an average thickness of 1.6 km, for a total of 2.85 million km3 of ice. While the ice sheet has waxed and waned during ice ages and warm periods, Greenland appears to have remained partially glaciated for at least 1 million years. Over the last century, the GIS has been volumetrically stable, with precipitation (adding volume) offset by melting and iceberg discharge (reducing volume). However, the GIS during the last two decades has lost about 280 km3 annually, which is the equivalent of 0.7 mm of sea level rise equivalent (SLRe) per year. The 2013 IPCC report (1) reviewed the evidence on the GIS. It concluded that several stable states of the GIS might exist; that the ice sheet might irreversibly shrink to a stable smaller state once a warming threshold is crossed for a certain amount of time; that the critical duration would depend on how far the temperature threshold has been exceeded and for how long; and that an irreversible decrease of the GIS appears very unlikely in the 21st century but is likely on multicentennial to millennial timescales in the largest warming scenarios. See SI Appendix, part C, for a further discussion. Modeling studies find a threshold temperature for GIS disintegration variously between 1 and 4 °C above mid–20th-century levels. However, it is misleading to suggest that complete GIS disintegration is inevitable when the temperature threshold is passed because the disintegration is relatively slow. Rather, as the IPCC (1) notes, “The complete loss of the ice sheet is not inevitable because it has a long time scale (tens of millennia near the threshold and a millennium or more for temperatures a few degrees above the threshold). If the surrounding temperatures decline before the ice sheet is eliminated, the ice sheet might regrow” (ref. 1, p. 1170) In thinking about tipping points for the GIS, it would be more accurate (although still oversimplified) to consider a threshold as measured in degree-years rather than degrees. The first question involves the equilibrium temperature–volume relationship. The paleoclimatic history of the GIS was thoroughly reviewed in Alley et al. (11). They summarize as follows: Paleoclimatic records show that the Greenland Ice Sheet consistently has lost mass in response to warming, and grown in response to cooling. [M]ajor changes of central regions of the ice sheet are thought to require centuries to millennia. The paleoclimatic record does not yet strongly constrain how rapidly a major shrinkage or nearly complete loss of the ice sheet could occur. The evidence suggests nearly total ice-sheet loss may result from warming of more than a few degrees above mean 20th-century values, but this threshold is poorly defined (perhaps as little as 2 °C or more than 7 °C [in regional temperature]). In their summary of the paleoclimatic record, Alley et al. (11) provide an equilibrium relationship between SLRe and global temperature, shown in Fig. 2. The current ice sheet has an equilibrium volume that, if completely melted, would add 7.2 m of SLRe. The paleoclimatic record has poor resolution of the global mean temperature at which the GIS will be nearly completely melted, but the interpretation of Alley et al. is that virtually complete melting will occur at between 1.3 and 4.5 °C sustained global warming. Other estimates show a similar level of warming and great uncertainty. For example, Robinson et al. (3) found a threshold for irreversible loss of 0.8–3.2 °C (this representing 95% confidence). The IPCC 1.5 °C report (2) concludes that the “evidence suggests that the temperature range of 1.5–2 °C may be regarded as representing moderate risk, in that it may trigger MISI [marine ice sheet instability] in Antarctica or irreversible loss of the Greenland ice sheet and it may be associated with sea level rise by as much as 1–2 m over a period of two centuries.” Fig. 2. Estimated equilibrium relationship between temperature change and equivalent sea level rise (SLRe) for the GIS. Data from ref. 11. The original figure is annual temperature over the GIS, and this is converted to global annual temperature using a ratio of 1.5:1 for GIS to global. The four points from Left to Right are modern (1960), the Holocene Thermal Maximum (HTM), the Eemian, and full ice sheet loss of MIS 11 or earlier, while the arrow suggests the impact of the glacial maximum, as discussed by Alley et al. This source also provides error bars on the estimates for the upper three points (from Left to Right) for temperature ±0.7, ±2.4, and ±2.5 °C and for sea level rise of ±0.5, ±2.0, and ±0.7 m. A key question is whether there are multiple locally stable temperature–volume equilibria for the GIS. The modeling studies of Ridley et al. (4) and Robinson et al. (3) find multiple equilibria for global temperature increases in the range of 0–2 °C. For example, Ridley et al. (4) simulate the long-run dynamics of the GIS with preindustrial forcings and different starting points from 0 to 100% of current volume. They find that the ice sheet volumes converge toward three stable equilibria at about 100%, 80%, and 20% of present-day volume (V0). As is common with dynamic systems, there are unstable equilibria between the stable ones; one unstable volume in Ridley et al. is around 90% of V0, which is the divide between the 80% and 100% of V0. The other unstable volume is around 50% of V0, which is the divide between 20% and 80% of V0. The study by Robinson et al. (3) examined the stability properties of the GIS for different temperature trajectories and found strikingly similar results to those of Ridley et al. This study finds hysteretic dynamics with multiple equilibria as shown in SI Appendix, Fig. G-1 and its legend. More precisely, they find three stable (and two unstable) equilibrium combinations of temperatures and volumes in the range of 0.6 to 1.5 °C. SI Appendix, Fig. G-2 shows the trajectory for ice-sheet volumes with different starting volumes and displays the hysteresis. According to ref. 3, at a global temperature increase of 1 °C, the separation point between the upper and middle equilibria is 80% of current volume, while the separation point between the middle and lower equilibria is about 40% of current volume. At low temperatures (<1/2 °C), all calculated paths go to the upper equilibrium volume, while at high temperatures (>2 °C), all paths go to the lower equilibrium volume. The GIS is therefore an example of a tipping system that might justify a temperature ceiling for policy. However, this point is where the integrated analysis of economics and geosciences becomes essential. Disintegration does not inevitably occur once the temperature threshold is passed. Rather, rapid and near-irreversible disintegration occurs only once a volumetric threshold is passed. A high-temperature path might well reduce the size of the GIS over, say, the next two centuries. However, as long as the GIS volume is above the volumetric threshold (say 80% of current volume), then reducing temperature back below its threshold will avoid passing the tipping point and prevent the resulting catastrophic melting.

Comparative Results of Alternative Ice Sheet Dynamics The paleoclimatic data do not yet provide a clear record for estimating the transient response of the GIS to different temperature trajectories. Understanding dynamics therefore relies on ice sheet modeling. As background, I examined the paths in several studies of GIS dynamics. Most studies take a trajectory for global or GIS warming and then track the ice sheet volume. A convenient way of summarizing the results is the melt rate per unit time per unit warming (in centimeters per century per degree Celsius). The results of the comparison are shown in Fig. 3. (For this discussion, the term “melt rate” is shorthand for the rate of decline of GIS volume per unit time.) Note that current models give highly divergent estimates of the transient response to warming. For example, Bindschadler et al. (12) conducted simulations for 500 y with seven ISMs and found a 500-y mean SLRe of 72.6 cm with a range of 8.5–142.6 cm. Also see refs. 13 and 14 and SI Appendix, Table C-2. Fig. 3. Alternative estimates of initial melt rate from different studies. Estimates are provided in SI Appendix, parts A and C. The arrow shows the range of studies for the Bindschadler et al. (12) model-comparison study for the 500-y horizon. [Legend: “Apple,” Applegate et al. (13); “Bind,” Bindschadler et al. (12); “Furst,” Furst et al. (14); “Ridley,” Ridley et al. (4); “Rob,” Robinson et al. (3). “DICE” is the results of the GIS estimates in the DICE-GIS model.] For calibration of the GIS dynamics, the DICE-GIS relies on the results from ref. 3, which have detailed simulations of GIS dynamics for different temperature trajectories. The numerical results of ref. 3 were provided by the authors, and projections for four of the temperature paths are shown in Fig. 4. The advantage of relying on this simulation is that the numerical results are available, and it has a wide range of temperatures as well as a long simulation period. The DICE-GIS calibration matches (3) primarily at high temperatures, as is seen in Fig. 3 but is in the middle of the pack of other studies. Fig. 4. Volume of GIS under different global temperature regimes from calculations in Robinson et al. (3). Data provided by A. Robinson (Complutense University of Madrid, Madrid, Spain).

Modeling the GIS for Inclusion in Integrated Assessment Models General Considerations. ISMs are highly complex as they require not only representations of the surrounding air and ocean temperatures but also, in the complete form, a 3D model of the dynamics of the ice sheet. The studies shown in Fig. 3 link climate models with ISMs. Such models allow changes in climate simulated by the climate models to interact with the ISM through surface mass balance (SMB) feedbacks. The feedbacks include changes in surface albedo and elevation, circulation changes induced by topographical change, and changes caused by changes in freshwater runoff. Since the modeling here relies on ref. 3, that study’s methods will be briefly described. The study starts with global climate models, which produce near-surface temperature anomalies prescribed over the boundary ocean points near Greenland. A regional energy–moisture balance model (REMBO) then takes the boundary conditions as well as the outputs of the ISM to simulate daily temperature and precipitation as well as SMB, snowpack thickness, and albedo. The REMBO outputs are inputs to simulation code for polythermal ice sheets (SICOPOLIS) model, which is a widely used 3D shallow-ice–approximation ISM. The relevant outputs of REMBO are SMB and surface temperature, which are inputs to SICOPOLIS; changes in topography and ice sheet extent calculated by the ISM are the output of SICOPOLIS and inputs to REMBO. The climate and SMB are updated every 10 ice sheet-model-years to provide accurate surface forcings to the ice sheet. Note that because REMBO is coupled to SICOPOLIS, the approach explicitly captures elevation and albedo feedbacks in the climate/ice sheet system at relatively high resolution (20-km grid). It is important to understand how the albedo–altitude feedback leads to instability. Warming will reduce the elevation of the ice sheet, which will lead to higher temperatures at the top of the ice sheet. Additionally, a warmer ice sheet will have less snow cover, reducing the albedo and adding further heat. For example, snow has an assumed albedo of as high as 0.8, while ice-free land has an assumed albedo of 0.2. Therefore, while only 20% of solar radiation would be absorbed by a cold ice sheet covered with snow, 80% of radiation would be absorbed by ice-free land. It is easily seen how this feedback could lead to continued deglaciation. The offset to this factor, it turns out, is the positive association of temperature and precipitation, which can offset some or all of the albedo–elevation feedback. Modeling Details. In developing the DICE-GIS model, it is necessary to find a numerical structure that represents GIS behavior in a robust and parsimonious manner. Call this a “reduced-form model.” The model must be simple enough to include in a few equations, yet reliable enough to represent the larger models. For example, the standard SICOPOLIS model has thousands of equations and clearly cannot be run in an optimization model. The strategy in developing the DICE-GIS model is to incorporate a simplified representation of more complex GIS models. The following presents a structural dynamic model that allows for any of the three types of dynamics (reversible, hysteretic, and irreversible). It is small, can be calibrated to larger realistic models, and can be incorporated into the DICE model. Eq. 1 above provided the (V, T) equilibrium. The next question involves the dynamics of volume adjustment. The simplest relationship is a differential equation (discretized in practice) in which the volume adjusts as a function of actual and equilibrium temperature and actual volume: ∂ V ( t ) ∂ t = g [ T ( t ) , T * ( t ) , V ( t ) ] . [2]The present study focuses primarily on the completely reversible system, shown in Fig. 1A. It initially assumes that the equilibrium function is linear to simplify the analysis (the other specifications are analyzed later). Eq. 2 is estimated from Robinson’s simulations using the data shown in Fig. 3 (see SI Appendix, part B, for the statistical results). The final equations are as follows: T * ( t ) = 3.4 [ 1 − V ( t ) / 100 ] , [3] Δ V ( t ) Δ t = − 0.0053 sgn [ T ( t ) − T * ( t ) ] [ T ( t ) − T * ( t ) ] 2 [ V ( t − 1 ) / 100 ] 0.2 . [4]Here, V*(t) and T*(t) are equilibrium volumes and temperature, while T(t) and V(t) are actual values. Eq. 3 takes the paleoclimatic equilibrium shown in Fig. 2 above, converts sea level rise to volume change, and linearizes the relationship between the modern era and the interglacial period. Note that the coefficient (3.4 °C) is the difference between the temperature at which the ice sheet is fully melted and the base temperature. At full volume (V = 100%), the equilibrium temperature is 0 °C, while the GIS fully melts in equilibrium at 3.4 °C global base temperature. Eq. 4 is the melt rate equation. The first term is the coefficient determined from a regression analysis (SI Appendix, part B). The second term introduces the sign of the temperature difference. The temperature difference enters as a quadratic function. The last term, [ V ( t − 1 ) / 100 ] 0.2 , ensures that volume is positive. To take an example, at an initial volume of 100% and a global temperature of 6 °C, the ice sheet decline is 0.19% per 5 y, or 28 cm of SLRe per century. If the actual temperature is less than equilibrium, the ice sheet rebuilds. Table 1, as well as SI Appendix, part B, show a comparison of the Robinson et al. (3) model calculations from with those of the DICE-GIS for high warming (6.7 °C). The ice sheet volume in the DICE-GIS calculation declines more rapidly at the beginning, but the two models have similar long-run trajectories. The differences between the Robinson and DICE runs are small relative to the differences among ISMs shown in Fig. 3. Table 1. Results for GIS volume for DICE and Robinson calculations for different temperature trajectories and time periods

Model Structure of DICE-GIS The DICE-GIS model is a straightforward integration of the reduced-form GIS model discussed in the last section with the DICE integrated assessment model. The DICE model is the latest version of a series of models of the economics of global warming developed at Yale University by Nordhaus and coworkers (15⇓–17). The model views climate change in the standard framework of economic growth theory known as the Ramsey model. The DICE model modifies the standard growth model to include climate abatement investments, which are analogous to conventional capital investments. The model contains all elements from economic growth to emissions to concentrations to climate change to damages in a form that attempts to represent simplified best practice in each area. A few changes have been introduced into the standard DICE model to reflect the long time period. The full set of variables and equations is provided in SI Appendix, parts D and E, and the changes in the DICE model for the GIS runs are included in SI Appendix, part E. The new assumptions in the standard DICE module are the following. First, no negative emissions are allowed past 2200. If these are allowed, then the optimal solution is to run atmospheric carbon concentrations low enough that the GIS stays at current volume. Second, the rate of decarbonization is set at zero after 2200. Without this assumption, emissions go quickly to zero. Third, several parameters are set as constants after 2200 for computational stability. These include the savings rate and the rate of productivity growth. The runs are for 1,500–5,000 y depending on the simulation. The following list shows the scenarios used for the present study. i) Discounting. Because of the long time lags, disintegration has a small impact on policy under normal discounting. The simplest way to deal with this concern is to consider a range of discount rates. This approach is consistent with other studies that advocate low discounting to reflect major losses in the distant future, such as the Stern Review (18).

ii) Damages on SLR. A second assumption concerns the damages from sea-level rise. The present study takes the results of Diaz (19). This study finds that SLR of 0.8 m in 2100 has an impact of 1.5% of global output without adaptation and 0.18% of output with adaptation. This study takes the intermediate estimate of 1% of global output lost for each 1 m of SLR. The damage function is linear in SLR in light of findings from Diaz. This function implies that complete disintegration of the GIS would lead to ≈7% loss in global income each year. Note that if modeling takes the constrained volume approach in iv, b below, the GIS component of damages is omitted.

iii) Alternative melt rates. The standard melt rate has been discussed above. For sensitivity analyses, I assume a melt rate two times the calibrated level. This takes the melt rate beyond any of the estimates that have been included in the 2013 IPCC report (1) but is useful for analytical purposes.

iv) Economic calculations. There are two alternative methods of treating the economic impact of the disintegration of the GIS. (a) Damage function approach. The first is to modify the standard DICE damage function to include GIS damages as described above in ii. (b) Constrained volume approach. A second approach is to constrain the GIS volume to be above a given threshold and remove the GIS economic damage function. For example, GIS volume might be constrained above 90% of its original volume. The volumetric approach is useful if estimates of GIS damages are imprecise, if the damage-function approach is too uncertain, or if it is desired to avoid tipping points in the ice sheet dynamics. This list provides a large array of potential strategies for including the GIS in integrated assessment models. Another set of issues is the potential for irreversibility and hysteresis, whose dynamics are examined in later sections.

Results for the DICE-GIS Model Major Results. Begin with the results of the standard DICE model with the GIS module added. Fig. 5 shows the trajectory of GIS volume for three cases: an optimal climate policy; a baseline of no climate policy; and a baseline policy followed by geoengineering after 500 y. Fig. 5. The evolution of the GIS under optimal and baseline (no-policy) cases plus a geoengineering scenario. See text for discussion. The arrow is the range of model estimates from different ISMs for a high warming scenario from IPCC (1), p. 1191, table 13-8, and has comparable forcings as the DICE-GIS baseline run. The optimal policy stays above the upper instability for the GIS found in studies discussed in text, whereas the baseline declines below the upper threshold. The run labeled “Base-geo” shows the result for a no-policy run followed by a geoengineering experiment after year 500 (see Geoengineering to Limit Temperature for details). The results are straightforward. With standard damages, discounting, and melt rate, the baseline path calculates that the GIS melts gradually over the coming centuries. By contrast, the optimal path shows a slower melt, and the GIS remains above 80% of current volume and is safely above the upper threshold volume found in refs. 3 and 4, as discussed above for at least a millennium. The results of the optimal and baseline run are shown in Fig. 5. The arrow in Fig. 5 shows IPCC estimates of the impact of no-policy radiative forcings on the GIS at 500 y for an array of ISMs (1). The figure also shows the result of a geoengineering experiment; this run assumes a baseline run for 500 y and then that a geoengineering technology reduces the temperature over the GIS to 0 °C starting at year 500. The GIS model used here suggests that the GIS rebuilds after geoengineering, but at a very slow pace. This point is further analyzed in Geoengineering to Limit Temperature. It is useful to compare the SCC for three cases: standard DICE damages only, GIS damages only, and combined damages. The SCC is defined as the marginal damage along the optimized path, which will also be the price of emissions on those paths in a market context. For runs that have temperature or volumetric constraints, these constraints are interpreted to represent implicit damage functions, but normal damages are always included as well. Table 2 shows that the marginal damages associated with GIS disintegration are a small fraction of the total damages, comprising 0.4% of the total. Further calculations in SI Appendix, part J examine the impact of alternative discount rates on the SCC. The easiest way to implement alternatives is to assume different constant discount rates, here 0.1%, 1%, 2%, 3%, 4%, and 5% per year and at two different melt rates. Inclusion of damages from GIS melting has a small impact in all cases. For the lowest discount rates and the higher melt rate, the GIS adds at most 5% to the estimated SCC. Table 2. Social cost of carbon This finding is important for the question of whether current estimates of the SCC underestimate the “true” SCC because of omissions of major tipping points such as the GIS. While the GIS is but one of the potential omissions in current methods, the calculations suggest that omitting the impact of GIS damages on the SCC is somewhere between small and negligible. Volume Constraints. Estimates of the damage from sea level rise due to GIS melt are highly uncertain. The damage estimate used in the modeling assumes limited adaptation, whereas, based on the study by Diaz (19), high adaptation would produce about one-fifth of the damages. An alternative approach is to limit the decline in the volume of the GIS. A natural set of limits would use the thresholds that have been suggested by current research. Based on the earlier discussion, the calculations here constrain the GIS volume to remain above 10%, 50%, and 90% of current volume. SI Appendix, Table J-2 shows the results of the volumetric approach for different discount rates and the standard and 2× melt rate (18 cases). The calculations provide estimates (i) for the standard non-GIS damage function plus GIS damages and (ii) for volumetric constrains and standard non-GIS damages. Case (i) therefore monetizes the GIS damages, while (ii) uses a physical constraint rather than monetary GIS damages. For most cases, the volumetric constraint is not binding, and the SCC is slightly lower than the standard estimates in Table 2. The volume constraint is binding in only four cases: For the upper threshold (90% minimum) and with the two higher discount rates (DICE and 3% discounting). Perhaps the most illuminating case is for DICE discounting (∼4.5% per year) and the base melt rate (call this the “standard case”). Here, the SCC rises from $31 to $39 per ton CO 2 —a SCC premium of $8. The reason for the increased SCC is that, for the optimal run in the standard case, the GIS volume declines below the 90% level. To reduce the melt rate sufficiently to remain above the 90% volume threshold requires raising the carbon price from $31 to $39 per ton CO 2 . This result indicates that quantitative targets, such as GIS volumetric constraints or other constraints to avoid tipping points, might be useful supplements to damage functions in cost–benefit studies. It should be noted, however, that the SCC premium is extremely sensitive to the volume constraint and changes sharply if the volume constraint changes slightly. Geoengineering to Limit Temperature. A final set of experiments examines geoengineering, which reduces global and GIS temperature to reverse the GIS disintegration. These experiments assume that global and GIS temperatures are reduced to 0 °C after a given year. The geoengineering might occur through radiation management (putting particles in the atmosphere) or carbon reduction (say through carbon removal technologies). In looking at these geoengineering simulations, the striking result is that rebuilding the GIS is quantitatively different from disintegration. The asymmetry is seen in Fig. 5, where the disequilibrium dynamics are very different in a melt mode from a rebuild mode. The asymmetry can also be seen in several geoengineering studies with large ISMs, including the simulations in both Ridley et al. (4) and Robinson et al. (3). Robinson et al. (3) shows a scenario in which temperature is reduced to 0.4 °C starting from a volume of 20% of current level. In this simulation, the ice sheet rebuilds to 70% of volume after 50,000 y. This result represents an increase of 0.07 mm/y. The estimates in the Ridley et al. (4) calculations indicate a buildup of about 0.1 mm/y in the accumulation phase. Applegate and Keller (20) estimate a buildup of 0.4 mm/y at 0 °C temperature (their figure 4). The geoengineering experiments in DICE-GIS are roughly the same as the results from the three modeling studies. Consider a scenario with a temperature anomaly of 6 °C global for 300 y, which leads to ∼2 m of SLR at that time. If the temperature is reduced to zero in a geoengineering experiment, the GIS is estimated to rebuild by only 0.2 m after 1,000 y, or about 0.2 mm/y. The reasons for the strong asymmetry is straightforward. An ice sheet melts when the temperature is elevated, and the decumulation (melting plus glacial discharge) exceeds accumulation (precipitation). However, there is no “negative melting” in the cold phase. Rather, to build an ice sheet requires not just cold temperatures, but also precipitation. When temperature declines, precipitation tends to asymptote to low levels, and the ice sheet buildup is consequently extremely slow. This implies that there is a sharp asymmetry in the response of the ice sheet to positive and negative temperature shocks. The conclusions of the geoengineering simulations have important implications for climate policy. The results apply to mitigation and carbon removal as well as solar-radiation management. They suggest that the disintegration of the GIS is essentially irreversible on a relevant societal timescale. The GIS will rebuild when temperature is reduced, but the growth is so slow that, from a human perspective, disintegration should be considered irreversible.

Alternative Equilibrium Specifications: Nonlinear, Irreversible, and Hysteretic The basic model analyzed here uses a linear relationship between equilibrium volume and temperature. This section considers alternative specifications of the equilibrium relationship: nonlinear, irreversible, and hysteretic. Nonlinear Equilibrium Function. The first alternative is to assume that the equilibrium volume–temperature function is nonlinear as shown by the finding from paleoclimatic studies and summarized in Fig. 2. Fig. 2 shows V* is a concave function of T*. The nonlinear relationship can be represented by a function of the form T ∗ = 3.4 ( 1 − V ∗ / 100 ) 0.5 . To maintain the same temperature–volume trajectory as the linear function for the first two centuries, the melt-rate coefficient is adjusted upward by about 10%. Calculations indicate that the optimal and baseline paths are virtually identical for linear and nonlinear specifications with standard coefficients and discount rates. For example, the SCC associated with only GIS damages is $0.134 per ton CO 2 with the linear model and $0.133 per ton CO 2 with the nonlinear model. Differences appear with low and superlow discount rates (1% and 0.1% per year). At these low rates, the ice sheet melts slightly more slowly with the nonlinear specification than the linear specification. This leads to slightly higher long-run volumes and lower SCCs with the nonlinear equilibrium function compared with the linear function. So the conclusion on introducing a nonlinear (concave) equilibrium temperature–volume relationship is that there are negligible changes for policy in the near term with standard parameters, while long-run disintegration is slightly lower with the nonlinear function. Numerical results are not presented as they are not interesting. Irreversible Disintegration. A second alternative structure assumes irreversible disintegration, as in Fig. 1C. That is, once melted to a given volume, the ice sheet cannot rebuild. To begin with, this is both physically and historically unrealistic. Paleoclimatic data (such as reported in Fig. 2) indicate a reversible pattern during ice ages and interglacial periods. Moreover, all models that allow for wide variations of forcings indicate changes from virtually ice free to highly glaciated conditions of the GIS. On the other hand, as the experiments with geoengineering indicate, the rebuilding of the GIS in a period of colder conditions is extremely slow. It is useful to determine, therefore, whether assuming irreversibility changes the DICE-GIS results. To test this question, a constraint was added that the change in ice sheet volume is always nonpositive for all alternative assumptions in the linear model. These runs show imposing irreversibility has no impact on the optimal or baseline policies that for all discount rates, two alternative melt rates, and the three target volume constraints. All these results indicate that, for the modeling of the ice sheet used in the present study, the dynamics of GIS melting is close to irreversible. The pace of rebuilding is so slow that, from a societal vantage point, the ice sheet dynamics can usefully be thought of as irreversible. Hysteresis in Equilibrium Temperature–Volume Relationship. A final approach is to examine the implications of an ice sheet displaying hysteresis. To modify the model for this property, the equilibrium temperature–volume relationship was assumed to follow a cubic function with the shape similar to Fig. 1B (the exact shape is shown in SI Appendix, Fig. H-1). The hysteretic curve is generated to resemble the estimates in Robinson et al. (3) shown in SI Appendix, Fig. G-2. The equation has an upper-branch tipping point at a temperature of 2.5 °C and a volume of 70% of current volume. The lower branch has a tipping point at 1 °C and 25% of current volume. The equation has an equilibrium of zero volume at a temperature of 3.4 °C. This cubic equation then replaced the linear equilibrium temperature–volume Eq. 3 in the standard DICE-GIS model. The melt rate equation is reestimated to fit the Robinson data. The major difference between the linear and hysteretic model is that the latter has a lower melt rate at the initial volume. As a result, for the early years and for a given temperature path, the GIS volume is higher with the hysteretic specification than with the linear specification. This result holds in all of the 16 variants of discount rates and melt rates for the optimal and baseline runs. This finding is inevitable if the two functions are constrained to have the same values at the two ends (100% volume at 0 °C volume and 0% volume at 3.4 °C). For the optimal policy, the results of the hysteretic model are very close to the linear model. The SCC for standard DICE parameters is 0.02% lower in the hysteretic case than in the linear case. The terminal volume (1,500 y out) is higher in the hysteretic case: 90.8% of original volume in the hysteretic case versus 85.4% in the linear case for standard parameters; and 97.8% versus 96.8% (hysteretic versus linear) of original volume for superlow discounting and the high melt rate. So to a first approximation, introducing hysteresis when policy is optimized makes little difference because policy ensures that melting stays away from the tipping point. For the baseline or no-policy case, the results are more illuminating. The SCC is lower in all cases for the hysteretic case (because of the lower initial melt rate). However, the ice sheet melts much faster once the tipping point is passed. With no policy, therefore, complete disintegration occurs more rapidly with hysteresis. Additionally, a hysteresis calculation was made using the lower threshold of 1.3 °C from Robinson et al. (3). There was no significant difference in the optimal policy or trajectory between the two hysteresis models. The results are described in SI Appendix, part H. So an important finding for the GIS is that the introduction of hysteresis instead of a monotonic (T*, V*) function makes little difference to the optimal policy. The reason is that the optimal policy stays away from the hysteretic threshold. However, for policies that pass the tipping point because of weak climate policies, hysteresis may make the outcome worse more quickly.

Conclusions and Qualifications The present study developed an integrated economic-geophysics model to analyze the interaction between the GIS and economic and energy factors. Integrating the DICE model with a small model of the GIS allows an internally consistent set of assumptions and results regarding emissions, climate change, GIS disintegration, sea level rise, and damages. While the different modules are simplified relative to high-resolution models, they have the advantage of integrating the different systems so that alternative policies can be assessed. The major message is that integrated modeling can improve our understanding of the complicated dynamic interaction of the economy and large-scale geophysical events as well as design policies to prevent crossing dangerous tipping points or irreversibilities. The major results were provided in the introductory section. This section concludes with a discussion of some of the qualifications with the current analysis, focusing on the major issues that arise from adding the ice sheet modeling. It leaves to the side standard issues of integrated assessment models such as DICE, which have been subject to extensive analyses. One major concern about including the GIS is that the equilibrium behavior is imperfectly understood. In particular, the answer to whether there are single or multiple equilibria is unclear. The best current evidence is that there are multiple equilibria for global temperatures increases between 0 and 4 °C. The evidence seems clear that virtually complete disintegration will eventually occur at global temperature increases above 6 °C, although “eventually” is many centuries. A second issue is the transitory dynamics of disintegration. As the survey above indicates, current models provide highly divergent estimates of the melt rate at different temperatures. A multimodel survey gives a range of estimates of melt rates of a factor of 7. The divergence arises because of the complexity of ice sheet dynamics and the absence of a precise transition history in the paleoclimatic record. A third uncertainty is the economic impact of sea level rise. The most careful study to date (19) indicates that there is a range of a factor of 10 for estimated impacts between a full-adaptation and a no-adaptation scenario. This uncertainty can be avoided by employing volumetric policies that constrain the disintegration of the GIS, but quantitative limits have the disadvantage of having a weak economic basis. A final issue that arises in analyzing the role of the giant ice sheets is the presence of large uncertainties (Fig. 3). An important question is how uncertainty would affect the outcomes and optimal policies. Preliminary estimates indicate that uncertainties about the melt rate and GIS damage coefficient have little impact on the SCC. Doing a full analysis of uncertainty is beyond the scope of the present study.

Acknowledgments The author is grateful for helpful comments from many researchers and colleagues. Particular thanks go to Klaus Keller, who has pioneered work on the interface of economics and earth sciences. Additionally, I am grateful for feedback from Richard Alley, Eli Fenichel, Ken Gillingham, Matt Kotchen, Robert Mendelsohn, Joe Shapiro, Martin Weitzman, as well as the editor and anonymous referees. The research was supported by the National Science Foundation and the US Department of Energy as well as by a fellowship from the Carnegie Foundation. All views and errors are the responsibility of the author.

Footnotes Author contributions: W.N. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

The author declares no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: A full set of files, sources, computer programs, and methods have been deposited in the Harvard Dataverse (https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/UYCPFL).

See Commentary on page 12134.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1814990116/-/DCSupplemental.