Modeling the Greenland Ice Sheet


  • Learn how to set up a coarse continental-scale Greenland model
  • Follow an example to initialize a continental domain, with a given ARGUS (*.exp) file and to parameterize with the SeaRISE netcdf dataset
  • Become familiar with how to set up and force transient input in ISSM
  • Plot results of forward simulation experiments

Go to trunk/examples/Greenland/ to do this tutorial.


In this tutorial, you will learn how to set up a continental Greenland model using the SeaRISE ice sheet model input dataset [Nowicki2013a]. In addition, you will gain experience in interpolation of datasets on to your continental ice sheet mesh and in setting up a transient forcing in ISSM. Finally, you will run a transient solution, resulting in a forward historical simulation of the Greenland Ice Sheet. Note that the model we set up here is coarse and is not recommended for use in a publication. A good use for this example it is use it as a starting point to learn how to use ISSM. You may wish to improve the model provided here by increasing the resolution of the ice sheet domain outline, increasing the mesh resolution, and choosing your own/improved datasets for model parameterization.

Tutorial steps to be taken:

  • Mesh Greenland with given *.exp file
  • Adapt mesh using SeaRISE velocity data
  • Parameterize (similar to the PIG model), except that all domain boundaries are on the ice front and do not have to be constrained
  • Stress Balance: run inverse method to control drag
  • Transient: Run a 20-year forward run
    • Use an appropriate time step for your resolution
    • Force SeaRISE surface mass balance for 10 years
    • For the next 10 years, simulate a warming scenario: decrease the surface mass balance linearly, reaching a decrease of 1.0 m/y by year 20
  • Plot transient results
  • Run an example exercise, forcing your Greenland model with historical SMB through time


In Step 1, we create a mesh using the triangle method (lines 10-11). This creates a new model named md and meshes the model domain, defined by an outline file 'DomainOutline.exp', at a resolution of 20,000 meters. Next, we adapt the mesh based on SeaRISE velocities, where the minimum resolution will be 5 km in locations where the velocity gradient is large and 40 km where the velocity gradient is small. The velicity data we will use resides in '../Data/' (line 5). Step 1 consists of the following steps:

  • Fill the variable vel with the interpolated velocities (Hint: you need x and y velocities plus ncfile x and y coordinates)
  • Mesh adapt your Greenland model (bamg)
    • Use variable vel
    • Set hmax=400000 and hmin=5000
  • Convert x, y coordinates to lat/long and then save your model to a file

Review the code used to create a continental Greenland mesh (lines 8-30) in the readme.m file. After creating the mesh and saving the model, the code uses plotmodel to plot a mesh visualization.

Execute step 1 in the runme.m file. After doing so, you should see the figure below:



Call the setmask function with empty arguments, to denote that all ice is grounded. Then parameterize your mesh with file Greenland.par. Next, set your flow equation to SSA for all. Read through the parameter file ./Greenland.par, which is similar to your PIG .par file, but for Greenland. Here, we are parameterizing a full continental domain, so all points along the domain boundary will be considered ice front. As a result, these boundaries do not need to be constrained, therefore the single point constraints variables will all be set to NaN.

Run step 2. This will save your parameterized model. Now, plot the new model thickness and velocity. For example:

>> plotmodel(md,'data',md.geometry.thickness) >> plotmodel(md,'data',md.initialization.vel,'caxis',[1e-1 1e4],'log',10)

Stress Balance

Use control methods to inversely solve for Greenland FrictionCoefficient (Step 3, lines 44-81). Note: Remember that md.inversion can be called for help!

  • Set three different cost functions
    • Absolute value of surface velocity
    • Log of surface velocity
    • Drag coefficient gradient
  • Set cost functions coefficients to 350, 0.6, and 2*10^-6
  • Set gradient scaling to 50
  • Specify max inversion parameter = 200, min inversion parameter = 1
  • Solve a 30-step Stress Balance model in 2D, SSA
  • Copy result Friction Coefficient to model (md) value
  • Save your model

Review step 3 in the runme.m to verify that the parameters have been set properly. Run step 3 in the runme.m to perform the steps above.


You are now ready to run a transient! In Step 4, we will simulate a simple constant warming trend over Greeland by forcing a temporal decrease in md.smb.mass_balance.


Specify a transient forcing by adding a time value to the end (in the end+1 position) of the column of forcing variable values. For example, let SMB be the initial values of surface mass balance. To impose the forcing such that before time 10, surface mass balance is set to the column vector smb, and after time 20, it is set to smb-1 we use the following code:

  • >> md.smb.mass_balance = [ smb smb-1] >> md.smb.mass_balance = [ md.smb.mass_balance; ... [10 20]]

By default, ISSM will linearly interpolate surface mass balance between time 10 and time 20 in this example. Prior to first and after last imposed time, forcing values remain constant. In order to turn interpolation off (i.e. use a step function), you would set md.timestepping.interp_forcings=0. If this values is set to 0, then your surface mass balance will change at the specified time, and will remain constant until a new value (column vector with time in the last row) is specified.

Steps to set up your transient:

  • Set control md.inversion.iscontrol back to 0
  • Interpolate surface mass balance from SeaRISE dataset, converting from water to ice equivalent
  • Impose SeaRISE surface mass balance for 10 years then linearly decrease to 1 m/yr by year 20
  • Set time step to 0.2 and output frequency to 1 (every time step will be output in results)
  • Ask your model to save IceVolume, TotalSmb, and SmbMassBalance transient output
  • Solve a 20 year Transient in 2D, SSA
  • Save your model
  • Review lines 83-112 in runme.m

In Step 5, we give you an example of how to plot the transient results (lines 114-145). To see how the transient results are stored in your model, type md.results.TransientSolution.

Now, run steps 4 and 5 to launch your transient and plot results.



Now, let's run our transient with historical mass balance! Use Jason Box's surface mass balance (SMB) time series as forcing [Box2013a,Box2013b,Box2013c].1

First, format the SMB provided. In Step 6 of the runme.m file, we extract the SMB timeseries from the netcdf file, and create a timeseries plot (lines 147-175). Execute step 6. This will result in the figure below:


In Step 7, we will relax the model towards equilibrium with the mean SMB forcing. An example of a 20 year relaxation to the time series mean is shown in runme.m, step 7. Run step 7, which will assign the mean SMB to md.smb.mass_balance and run a transient model for 20 years, with a timestep of 0.2 years, saving the results every timestep. Step 7 will save the results in the Model "Greenland.HistoricTransient."

To plot the relaxed version of the model that you just created, change step 5 to load the model "Greenland.HistoricTransient" rather than "Greenland.Transient," and run step 5 again.


To reach equilibrium, the model should run on the order of 1000 years. Since, 1000 years might take quite a long time to run on a personal computer, you may want to try running for 200 years instead.

To accomplish this extended relaxation, alter step 7 to run for the extended time period (200 years instead of 20 years). In the last line of this step, save your model as ./Models/Greenland.HistoricTransient_200yr instead of ./Models/Greenland.HistoricTransient, to avoid overwriting the old model. Then, run step 7 again. This run of 200 years will take longer than your orginal 20 year run.

When you are done with step 7, complete step 8 on your own as an exercise. Fill in the required code to plot the results in step 8. Follow the comments, write the code to load the historic transient model, and create line plots of relaxation run (use Step 5 as a reference). Then, save surface mass balance by looping through 200 years (i.e. 1000 steps). Plot the surface mass balance time series in the first subplot. Title this plot "Mean Surface Mass Balance".

Next, save velocity by looping through 1000 steps. Plot velocity time series in a second subplot. Title this plot "Mean Velocity".

Lastly, save Ice Volume by looping through 1000 steps. Plot volume time series in a third subplot. Title this plot "Ice Volume" and add an x label of "years". The resulting plot should look like this:


In Step 9, we will use the 200 year relaxed ice sheet as a starting condition for a historic transient run. To do so we need to save the 200 year resulting geometry and velocities into the model state. To load your past results see lines 254-259.

Next, we load the Box time series saved earlier in mat file, smbbox.mat, and then (lines 261-300):

  • Interpolate every month of Box SMB onto the ISSM grid: insert a column for each month
  • Add a final row indicating that the value should be set in the middle of each month
  • Solve at a monthly time step and save monthly results

Run step 9, which will excute your historical transient forward simulation, monthly from 2003-2012.

Then, run step 10 to plot a time series of total surface mass balance, max velocity, and ice volume. See Lines 305-329.


1 The year 1840-2012 Greenland near surface air temperature (T) and land ice SMB reconstruction after Box [2013] is calibrated to RACMO2 output [Meijgaard2008,Ettema2009,Broeke2009,Angelen2011]. The calibration for T and SMB components is based on the 53 year overlap period 1960-2012. The calibration for snow accumulation rate is shorter because ice core data availability drops after 1999. Calibration is made using linear regression coefficients for 5 km grid cells that match the average of the reconstruction to RACMO2.

The RACMO2 data are resampled and reprojected from the native 0.1 deg (LaTeX equation10 km) grid to a 5 km grid better resolving areas where sharp gradients occur, especially near the ice margin where mass fluxes are largest. Several refinements are made to the Box [2013] temperature (T) and SMB reconstruction. Multiple station records now contribute to the near surface air temperature for each given year, month and grid cell in the domain while in Box [2013], data from the single highest correlating station yielded the reconstructed value. The estimation of values is made for a domain that includes land, sea, and ice. Box [2013] reconstructed T over only ice. A physically-based meltwater retention scheme [Pfeffer1990,Pfeffer1991] replaces the simpler approach used by Box [2013]. The RACMO2 data have a higher native resolution of 11 km as compared to the 24 km Polar MM5 data used by Box [2013] for air temperatures. The revised surface mass balance data end two years later in year 2012. The annual accumulation rates from ice cores are dispersed into a monthly temporal resolution by weighting the monthly fraction of the annual total for each grid cell in the domain evaluated using a 1960-2012 RACMO2 data.

Additional Exercises

  • Increase SMB instead of decrease over time
  • Create an instantaneous step in SMB forcing at 10 years instead of a steady change over time
  • Create a more advanced SMB forcing, like cyclic steps or a curve
  • Force SMB to change only in certain areas of the ice sheet
  • Add more melt in the ablation zone, but more snow in the upper elevations
  • Force another field transiently (e.g. friction coefficient)
  • Run the Box time series yearly or for a longer subset of time. This could take a while!


  • J. E. Box, N. Cressie, D. H. Bromwich, J-H. Jung, M. van den Broeke, J. H. van Angelen, R. R. Forster, C. Miege, E. Mosley-Thompson, B. Vinther, and J. R. McConnell. Greenland Ice Sheet Mass Balance Reconstruction. Part I: Net Snow Accumulation (1600-2009). J. Clim., 26(11):3919--3934, JUN 2013.
  • Jason E. Box. Greenland Ice Sheet Mass Balance Reconstruction. Part II: Surface Mass Balance (1840-2010). J. Clim., 26(18):6974--6989, SEP 2013.
  • Jason E. Box and William Colgan. Greenland Ice Sheet Mass Balance Reconstruction. Part III: Marine Ice Loss and Total Mass Balance (1840-2010). J. Clim., 26(18):6990--7002, SEP 2013.
  • J. Ettema, M. R. van den Broeke, E. van Meijgaard, W. J. van de Berg, J. L. Bamber, J. E. Box, and R. C. Bales. Higher surface mass balance of the Greenland Ice Sheet revealed by high-resolution climate modeling. Geophys. Res. Lett., 36:1--5, JUN 16 2009.
  • S. Nowicki, R.A. Bindschadler, A. Abe-Ouchi, A. Aschwanden, E. Bueler, H. Choi, J. Fastook, G. Granzow, R. Greve, G. Gutowski, U. Herzfeld, C. Jackson, J. Johnson, C. Khroulev, E. Larour, A. Levermann, W.H. Lipscomb, M.A. Martin, M. Morlighem, B.R. Parizek, D. Pollard, S.F. Price, D. Ren, E. Rignot, F. Saito, T. Sato, H. Seddik, H. Seroussi, K. Takahashi, R. Walker, and W.L. Wang. Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project II: Greenland. J. Geophys. Res., 118:1--20, 2013.
  • W.T. Pfeffer, T.H. Illangasekare, and M.F. Meier. Analysis and modeling of melt-water refreezing in dry snow. J. Glaciol., 36(123):238--246, 1990.
  • W.T. Pfeffer, M.F. Meier, and T.H. Illangasekare. Retention of Greenland Runoff by Refreezing: Implications for Projected Future Sea-Level Rise. J. Geophys. Res. - Oceans, 96(C12):22117--22124, DEC 15 1991.
  • J. H. van Angelen, M. R. van den Broeke, and W. J. van de Berg. Momentum budget of the atmospheric boundary layer over the Greenland ice sheet and its surrounding seas. J. Geophys. Res. - Atmospheres, 116:1--14, MAY 18 2011.
  • Michiel van den Broeke, Jonathan Bamber, Janneke Ettema, Eric Rignot, Ernst Schrama, Willem Jan van de Berg, Erik van Meijgaard, Isabella Velicogna, and Bert Wouters. Partitioning Recent Greenland Mass Loss. Science, 326(5955):984--986, NOV 13 2009.
  • E. van Meijgaard, L. H. van Ulft, W. J. Van de Berg, F. C. Bosvelt, B. J. J. M. Van den Hurk, G. Lenderink, and A. P. Siebesma. The knmi regional atmospheric model racmo version 2.1, technical report 302. Technical report, KNMI, De Bilt, The Netherlands, 2008.