Simulating an ice sheet

Most of our future projections about the effects of climate change are informed by computer models. These models are based on the laws of physics and aim to capture a part of the climate system (e.g. the atmosphere, the ocean, the ice sheets). They can stand alone, or they can be coupled to other models to better represent the feedback loops that often exist at the boundaries of Earth's systems. Atmospheric warming results in a warmer ocean, which likely affects ocean circulation (such as the thermohaline pump) and currents (like El Niño and La Niña), all of which accelerate the melting of ice sheets, which reduces the amount of sunlight reflected, in turn warming the atmosphere more... you get the point. It's all interconnected.

In this post we'll be focusing on a single part of this system: the Greenland ice sheet. It's the second largest body of ice in the world next to Antarctica, but it's melting the fastest: it's estimated to have lost around $243\text{Gt}$ per year over the period 2010-2019. In contrast, the Antarctic sheet's estimated loss rate is "only" $148\text{Gt}{\text{yr}}^{-1}$ over that same period.¹

We'll be using the SICOPOLIS model to simulate the evolution of the Greenland ice-sheet under a particular high-emissions scenario. But what is a model, and why should you care?

A model is a simplification of reality. When we neglect air resistance, assume a spherical cow, or say that a car uses a certain amount of fuel per distance traveled (even though this depends on the speed of the car), we are creating a model. This means that they're fundamentally ignoring things by not taking them into account, but as the famous quote attributed to George Box goes:

All models are wrong, but some are useful.

The thing that makes models "wrong" is also what makes them useful most of the time. It would take centuries to simulate every single atom in the Greenland sheet, but realizing that the forces that dominate glacier movements happen on a much larger scale is the key that allows us to simplify. Of course, these are assumptions that you should test, but that's part of the scientific process. One way of validating models is by having them "predict the past", which we can compare to historical data like the geological record.

Most computer models operate using a fundamental equation:

$$\begin{array}{lcr}& \overrightarrow{𝐱}(t_{i+1})=ℒ(\overrightarrow{𝐱}(t_i))& \text{<span class="tml-eqn"></span>}\end{array}$$

This is saying that our state $\overrightarrow{𝐱}$ at the next timestep $t_{i+1}$ is given by some function $ℒ$ (our model) applied to the state at the current timestep $t_i$. In ice-sheet models, that function $ℒ$ contains 2 main physical processes (usually in the form of partial differential equations):

1. Snow fall and melt (which together make up the surface mass balance)
2. Ice flow; the physical movement of the glacier (internal deformation, basal sliding, calving)

Some models will include more or less detail based on the assumptions they make (mostly related to internal stresses): the shallow ice (SIA) and shallow shelf (SSA) approximation being the most common. I won't go into more depth on them here but if you want to read more about these processes and assumptions, I recommend the antarcticglaciers.org page: they have great sources on glacier flow and mass balance, for example. You might also try this great open textbook.

Apart from these processes, models also rely on external inputs. This can be temperature, precipitation or sea level data, but also the initial state of the ice sheet. These are the knobs we can turn, after which we let the physics play it out. And it's with this step in mind that we can start modeling our ice sheet!

We're going to start by installing SICOPOLIS. Installing NetCDF can be quite a hassle on MacOS (which I am); I recommend using MacPorts in that case. Note that we won't really need LIS (we're only simulating land ice), so you can skip that part of the installation if you want.

After it's installed, we're going to retrieve the input data for our model run. We're going to need to tell the model what the temperature is going to be in the future. I'm going to use one of the IPCC's Representative Concentration Pathways (RCP), which is a scenario of our future CO₂ emissions. I'm going to go with the quite pessimistic emissions scenario RCP 8.5, so we can really see the difference in the ice sheet. The original RCP scenarios are only defined until 2100, so we're going to use ECP 8.5 (Extended Concentration Pathway), which extends until 2300.²

RCP and ECP emissions scenarios. Reproduced from 2.

This is, however, only emissions data: we don't know what effect this will have on temperature. Luckily, people have taken these emissions scenarios and run them through computer simulations for us. Thanks to CMIP5, these simulations produce standardized data that we get to compare and browse in a single place. I personally prefer the Copernicus Climate Data Store, so head over there and have a look at the CMIP5 monthly data on single levels dataset.

Under the "Download data" tab, we'll pick RCP8.5 (the scenario we chose), "2m temperature" as our only variable (if you want to take this simulation a step further, you could also take the precipitation data and feed it into SICOPOLIS), and for the "Model" you can pick any (I went with IPSL-CM5A-LR because its projections are relatively close to the median). The other fields should either auto-select or have only a single option, which you should pick. Finally, click the "Show API request" button to generate some code with your request, which we'll paste into our editor and modify like this:

py

Before this will work though, follow the instructions here to setup cdsapi and accept the Terms of Use for the dataset. After doing this, your request should go through, producing a download.zip file which, after unzipping:

py

...gives us a netCDF file that looks a little like this.

py

We now have a file containing gridded temperature data for the whole Earth from the year 2006 to 2300. We could just globally average this, but to be a little more accurate I'll be averaging only the data that contains Greenland. To do this we'll make a mask for data points between latitudes 59° and 83°N and longitudes 11° and 74°W. In the lat and lon format of our data, that's bounds of lat: [59, 83] and lon: [286, 349]. To ignore points that aren't land, I'll filter out temperatures warmer than 0° C in the final time-slice of our data (we can assume Greenland won't move much in 300 years so this mask should be valid whenever).

py

Averaging these temperatures over time then gives a curve that we'll smooth a little bit using a rolling mean:

py

Nice! We now have the temperature increase with respect to the starting year 2006. What's left is to offset this curve by the actual temperature it was in 2006, and then we write this out to a file SICOPOLIS understands. To do this, have a look at the file sico_in/general/grl_warming_scenario_rcp26.asc. This file contains Greenland surface temperature anomalies (which is strictly speaking also what we've been dealing with: departures from a reference temperature, which we've set at 0). The file format is a tad undocumented, but the first line contains a comment specifying # <START_YEAR_OFFSET> <TIMESTEP> <END_YEAR_OFFSET>, all times and values being relative to 1990 (the common start year for SICOPOLIS runs). Each line then contains a year and a temperature value: the first line has the year -140, which relative to 1990 is 1850, and the comment at the bottom also confirms this.

Our current temperature anomalies are relative to the year 2006, but we should make it relative to 1990. To do this, we'll simply find the anomaly during 2006 in the file and add it to our temperatures. We'll also make the years relative to 1990.

py

We'll write this out to a file so we can use it later.

py

Now that we have our inputs done, we can actually run the model. SICOPOLIS runs are specified by so-called "header" files, which you'll find in the headers directory. We're going to base our run on the existing file sico_specs_repo_grl10_b2_future21_asmb.h. I'm changing the resolution from $10$ to $20$ km to speed up our simulation (we can run the resolution doubler afterwards if necessary for plotting), and also the start and end time, to make sure it runs during the range of our data (2006 - 2300). I'm also setting GRIP_TEMP_FILE to our newly created temperature anomalies, and specifying the time-slices at which I want data to be written out using TIME_OUT0. The changes I made are these:

f90
sico_specs_repo_grl10_b2_future21_asmb.h
sico_specs_grl20_rcp85.h

Since our initial conditions depend on the result of another simulation (as specified by ANF_DAT=3), we'll have to run that one first. This script will run both models for you:

sh

Once that is done, your output files will be in the folder sico_out/grl20_rcp85 (or whatever name you gave to your specfile). An explanation of what these different files contain can be found here, but the gist of it is that the .ser or _ser.nc files contain scalar variables (ser meaning time series), the .site or _site.nc also contain (different) scalar variables sampled at particular sites (which are predefined and specific to your modeling domain), and the files that end in 0001.nc, 0002.nc etc. contain full 2D/3D fields (e.g. ice velocity or thickness). These fields are written out at the times we set on line 1322 of our specfile, so in our case it's once every 50 years.

Your directory should look something like this:

As a first sanity check, we'll get the surface temperature anomaly from the grl20_rcp85_ser.nc file to make sure that it corresponds to the input data we generated. I'll also calculate the expected sea level rise, because why not.

py

The temperature anomaly matches our input, so that worked. We see that under this ECP8.5 scenario we can expect a sea level rise of around 4 m by 2300, which for me as a Dutchman isn't great news. In order to create our final plots, I'll do some more processing. The code is a bit lengthy so I've hidden it to make the text flow better, but you can view it below if you want. It will also be visible in the notebook if you choose to run it yourself.

This content was hidden. Click to reveal it.

Show

py

py

py

py

Here we can see the surface mass balance (SMB) of the glacier over time. Remember how we talked about snow fall and melt? It turns out there's a lot more ways for a glacier to gain and lose mass, processes which we respectively call accumulation and ablation. When you add these together, you get the net annual change in mass. The equilibrium line (where SMB = 0) is the part where this surface mass will "stay the same" year-to-year. It separates the areas where there is a net gain (accumulation) and loss (ablation). After the year 2206, there's no accumulation areas left on the glacier: we have a net surface mass loss at every point. At some point between 2206 and 2301 the lower part of the ice sheet even detaches from the main body. As the ice sheet retreats, it also loses contact with the sea - this will be important later. Let's now convert these maps into something we can more easily visualize over time.

py

py

Here I've plotted three scalar quantities. The integrated SMB is exactly that: the total sum of the SMB at every point in our glacier. We see that in this climate scenario, the ice sheet is rapidly losing volume. We see a slight recovery around the year 2250 - this might be because in our smoothed temperature data, we inadvertently modeled a slight "overshoot" scenario at the end (it drops a little), or because the ice sheet has lost contact with the sea completely (removing calving and warming through seawater as drivers). These hypotheses would require more investigation, however. Note that the Greenland ice sheet has a response time of around 1500-3000 years, so this is nowhere near the steady state it's going to eventually reach.

The ice discharge (D), in this case, is the mass lost to the ocean through calving. As the glacier melts, its terminus retreats, increasing the distance to the ocean, and so the discharge decreases. It's important to note that this result is wildly incorrect: it's underestimating the observed discharge by at least a factor of 10.³ This could be because of our relatively large grid size of 20 km, because other similar simulations with smaller grid sizes do get this right, but I'm not sure. It's important to investigate the limitations of your models & where (if) your assumptions break down. If this was a scientific publication, I would either redo the entire simulation and ignore it ever happened, or push this work to other people by mentioning it in a cop-out "Future work & recommendations" section.

The accumulation-area ratio (AAR) is given by the area where accumulation takes place divided by the total glacier area. For non-calving glaciers, this is often between 0.5 and 0.6, but for calving glaciers it can be around 1. Though a useful proxy, it doesn't tell us everything about the health of a glacier: we can have an AAR of 1 (net snowfall everywhere), yet the mass loss through calving could be so large that the glacier is still shrinking in volume. The SMB contributes to a clearer picture, but the climatic mass balance⁴ (the sum of SMB and the internal mass balance IMB) tells us what the whole glacier volume is doing over time. The IMB contains a lot of different internal processes which are out of scope for this article (for now).

py

py

Finally, here are 3 more variables (plotted for the year 2106) that can give us some insight into the horizontal movement of the Greenland ice sheet. We see that most of the ice sheet moves very slowly, at rates maxing out at 100 m $\cdot$ yr${}^{-1}$, but some areas move at rates up to 7 times that! These areas also have a nonzero slip ratio (the ratio of basal velocity $v_b$ over surface velocity $v_s$). At these areas where $v_b/v_s>0$, we get basal sliding: the base of the ice sheet is no longer stuck to the ground, but sliding along it. This is facilitated by a warm ice base, which usually produces a layer of water that allows the ice to slide across the rocky terrain. We can see this in the bed temperature plot: cold areas don't slide, and the edges (where the ice approaches its melting point) do. These fast-moving areas tend to have large slopes, and often end up being or leading to so-called outlet glaciers where calving occurs.

If you managed to get to this point, thanks for reading! I hope you learned a little bit about the basics of glaciers processes & how we model them. If you're hungry for more, I can recommend this more in-depth study, which also mentions other software and their respective assumptions and differences. I hope we can also agree that this climate scenario is not really a desirable future, so let's like, not get there. Maybe take the bike instead of the car. It's also healthier.