Home

GRenoble Ice-Shelf and Land-Ice model: a practical user guide

1. 1 Snow if itis present is melted at rate Csnow Melt water is supposed to percolate into the snow layer and refreezes as superimposed ice Runoff occurs when the amount of superimposed ice exceeds 60 of the yearly snowfall 2 Superimposed ice is melted at rate Cice 3 Glacier ice is melted Depending on the available energy not all steps or even none of the steps will be carried out The melting factors C and Cice are typically set to 0 003 and 0 008 mm d C Csnow and Cice take different values to account for the difference in albedo between snow and ice The portion of melted water that can refreeze denoted by csi is usually set to 60 but can reach values up to 70 Pfeffer et al 1991 ISOSTATIC RESPONSE The bedrock reacts to changes in ice load Figure 8 In GRISLI this isostatic response is described by the ELRA elastic lithosphere relaxed asthenosphere method Le Meur and Huybrechts 1996 The lithosphere is treated as a thin elastic plate This allows a regional response to the ice load so that the bedrock in a certain distance from the imposed ice load can still be affected in contrast to the rather simple local lithosphere method which accounts only for a local response The underlying asthenosphere is assumed to be viscous Due to its viscous property the asthenosphere responds to an imposed ice load with a time lag the characteristic relaxation time 7 In GRISLI the default value of this parameter
2. Parametres du run Hemin40 amp runpar nom du bloc parametres du run runname shelf002 icompteur 0 reprcptr d431Apvi5 k0100 CPTR pour depart a 0 itracebug 0 num_tracebug 166 comment_run SPIN UP_1 NH40_035 200k SPIN UP_2 NH40_41b 200k runname nom de 1 experience 8 caracteres icompteur reprise dans un fichier 0 gt non 1 gt oui 2 gt T et Hwat f 3 gt T seulement reprcptr nom du fichier amp netcdf bloc netcdf dtncdf 10 dtncdf espace en temps entre les sorties amp grdline bloc grounding line igrdline 0 if GRenoble Ice Shelf and Land Ice model igrdline 1 ligne d echouage fix e sinon 0 f amp timesteps bloc timestep tend 30000 tbegin 0 dtmin 2 e 2 dtmax dos dtt Bis 15 testdiag 0 02 0 016 tous les temps en annees dtt pas de temps long tbegin et tend debut et fin du run pour equation masse pas de temps mini gt dtmin maxi gt dtmax testdiag pour g rer le pas de temps dynamique dt ordres de grandeur 40 km dtmin 2 e 2 a moduler selon dx dtmax 1 dtt 5 tesdiag 0 02 amp topo_file filel NH_40km_140k_topo2_withice_GRISLI xyz file2 NH_40km_140k_topo2_withice_GRISLI xyz filecoord Coordinates_NH_40km_GRISLI xyz filel topo de depart topo 21k g40 file2 topo de reference hemin2 g40 coordfile grid coordinates
3. m itracebug switch to debugging mode 0 gt no 1 gt yes num_tracebug corresponds to the number of the call within the routine for debugging The call to the debugging function is used to to print some matrix and arrays during run time and check the values and dimensions E comment_run character chain written in the ASCII output file and set by the user to write specific comment about the simulation Netcdf output frequency amp netcdf bloc netcdf dtncdf 10 dtncdf espac n temps entre les sorties This block is overwritten by the file GRISLI SOURCES Fichiers parametres hemin40 TEMPS NETCDF dat where the frequency of output is defined See Section 8 The user may therefore set an ad hoc value Grounded line physics amp grdline bloc grounding line igrdline 0 igrdline fixed gt 1 free gt 0 This blocks determines whether or not the grounded line will be fixed e not dynamic or will be able to move during the simulation Fixed grounded line is typically used to spin up Greenland to present day topography for example and avoid that the expansion of the ice sheet beyond the observed limit This has consequences on the ice shelves and on the mass conservation Model time steps Centro Euro Mediterraneo sui Cambiamenti Climatici SJ Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers amp timesteps bloc timestep tend 3000
4. w RS E N amp amp PS PH YE Geothermal heat flux mW m 2 Figure 4 Map of geothermal heat flux mW m2 from Shapiro and Ritzwoller 2004 interpolated on the Eurasian 20 km grid Climate correction during run time Similarly to the off line downscaling procedure changes in ice sheet elevation during run time affect input temperature and precipitation In GRISLI temperatures account for elevation changes at each GRenoble Ice Shelf and Land Ice model a practical user guide time step by means of a lapse rate value Teor t To A S t So 4 where Tor and To are the corrected and initial temperature fields respectively S is the elevation at time step t and Sp is the initial elevation The lapse rate A as described above takes different values for annual mean and July temperatures Precipitation is also corrected for changes in elevation Peor t Po exp y Teor t To 5 where Poo and Po stand for the corrected and initial precipitation respectively 6 GRISLI boundary conditions As boundary condition GRISLI uses a map of present day geothermal heat flux GHF recon structed by Shapiro and Ritzwoller 2004 Figure 4a The original GHF file is available at http ciei colorado edu nshapiro MODEL ASC_VERSIOW It has to be interpolated onto GRISLI grid before running and written in an ASCII file located in GRISLI INPUT HEMIN40 for example 1 72 5903244 2 1 138195953 3 L 76 1473999 241
5. h Centro Euro Mediterraneo sui Cambiamenti Climatici N Centro Euro Mediterraneo sui Cambiamenti Climatici N CMCC Research Papers o Catabatic Wind Precipitation ed Y a N ee e Heat Ice formation Sea ice Iceberg calving O CONTINENTAL SHELF Figure 6 Sketch illustrating the different regions considered in GRISLI ice sheets ice streams and ice shelves Figure from H Grobe AWI FUNDAMENTAL EQUATIONS The governing equations of ice sheet dynamics are based on the conservation of mass momentum and energy Op O ae V ou 0 6 u Pa V T p8 7 OT pc F u vr V KVT Qi 8 where p stands for the ice density u for the ice velocity 7 for the stress tensor g for the gravitational acceleration c for the heat capacity of ice T for the ice temperature for the thermal conductivity of ice and Q for the deformational heat If we consider ice as an incompressible fluid with a constant density the continuity equation 6 simplifies to vena 2 Ov Ow o 9 This equation is used to diagnose the vertical velocity from horizontal velocities the computation of the horizontal velocity will be described later on The vertical integration of equation 6 from the base of the ice sheet B to the surface elevation S leads to an equation for the ice thickness Defining the ice thickness H as surface elevation minus the base of the ice sheet S B assuming that the ic
6. 2 the basal melting occurring at a depth above the defined limit is a function of surface temperature set as a coefficient bmelt coe ftemp bmshel f 41 If the depth is larger then the defined value then bmelt bmshel f 42 Centro Euro Mediterraneo sui Cambiamenti Climatici CO Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers 3 At the grounded line for the point located along this line but which are still floating the basal melting is also set as a function of surface temperature since the depth of the grounded line might vary in function of the area bmelt coe ftemp bmgrz 43 If the point is at the grounded line and is not floating the basal melting is a combination between the grounded zone value bmgrz and the basal melting of the floating part bmelt ngr 4 x bmgrz x coe ftemp 1 ngr 4 x bmshel f x coe ftemp 44 where ngr corresponds to the number of grounded neighbour points In the last version of GRISLI no difference is made for the grounded zone and central ice shelf basal melting values if the grid domain is not Antarctica The user may easily edit the basal melting routine to account for different values between those two areas SURFACE MASS BALANCE The surface mass balance SMB in the ice sheet model is composed of accumulation and ablation SMB Accumulation Ablation The accumulation corresponds to the mean annual total precipitation obtained from
7. water level threshold to allow for ice stream over the sediment area to be combined with seuil_sedim GRenoble Ice Shelf and Land Ice model a practical user guide seuil_sedim minimum sediment thickness above which ice stream are triggered if seuil_hwater is reached concomitantly seuil_neff minimum effective pressure to trigger ice streams in Pascal Geothermal heat flux amp ghflux fileghf ijphi_hemin40 dat coefghf 1 geothermal heat fluxes fileghf input geothermal heat fluxes map coefghf to modulate the amplitud of the geothermal heat fluxes In this block the user can prescribed the geothermal heat flux file As for the sediment map it has to be interpolated on the GRISLI grid domain m fileghf input file of geothermal heat flux located under GRISLI INPUT HEMIN40 E coefghf coefficient set to modulate the geothermal heat fluxes coefghf 1 implies no changes in the input geothermal heat flux values In the routine ghf ghf x coefgfh Calving amp calving l nom du bloc calving m thode Vincent Hcoup 200 ifrange 3 meth_Hcoup 0 Hcoup epaisseur de coupure ifrange 0 gt NO CALVING ifrange 1 gt traitement de Vincent avec ice shelves frangeants biased ifrange 2 gt ice shelves frangeant seulement si bm bmelt positif biased ifrange 3 gt option 2 corrected but no test on ice thickness of upstream point ifra
8. Dim_hemin40 mod_communs mod_ell Liste_hemin40 diagnoshelf Liste_Netcdf routines_communes routine_elliptiques L NCDF_LIB lnetcdf L NCDF_LIB lnetcdff Liste_BLAS In the case of using the IBM Fortran and C compilers the BLAS libraries are included through the mk1 compilation option and the user has to modify the compilation block such that Hemin 40 Dim_hemin40 mod_ communs mod_ell S Liste_hemin40 diagnoshelf Liste_Netcdf routines_ communes Y routine_ elliptiques S LK o RUND Hemin 40 Dim_hemin40 mod_ communs mod_ell Liste_hemin40 diagnoshelf Liste_Netcdf routines communes routine _elliptiques L NCDF_LIB lnetcdf L NCDF_LIB lnetcdff Centro Euro Mediterraneo sui Cambiamenti Climatici D CMCC Research Papers Finally to compile and run the model for the Northern Hemisphere domain gt cd GRISLI SOURCES gt make Hemin 40 This will create and executable located in GRISLI bin gt cd bin gt Hemin 40 The other grids included with GRISLI are listed at the beginning of the Makefile 4 GRISLI grid domain The equations governing ice flow ice thickness and ice temperature evolution are discretised with the finite difference method The prognostic variables are computed on a rectangular regular grid The grids are created off line using any spherical or ellipsoidal geographical projection of the topography decided
9. edit the module choix The Makefile is structured as follows NetCDF libraries compilers and compiling options The user has to modify those options to adapt them to his own operative system and compiler NETCDF libraries NCDF_INC usr local netcdf 4 1 2_ifort include NCDF_LIB usr local netcdf 4 1 2_ifort lib Compiler IFORT ifort Compiling options and FLAGS ARITHM O2 xHost no prec div mkl m64 debug g CB fp stack check check all FT IFORT ARITHM c IS NCDF_INC g LK IFORT ARITHM 9 F_NETCDF S IFORT ARITHM c IS NCDF_ INC The list of the routines and modules that are compiled whatever is the grid domain mod_communs runparam mod o 3D physique gen_mod o pdd_declar_mod o iso_declar_mod 0 3 o isostasie mod 0 3 0 noisostasie_mod 0 3 0 routines_communes tracebug o ablation_Thompson_1997 0 ablation annual month o The routines associated with the selected grid domain When implemented a new grid the user has to duplicate this part of the Makefile to create the list of routine associated with his new grid domain GRenoble Ice Shelf and Land Ice model a practical user guide Liste_hemin40 output_hemin40_mod 0 5 0 lect clim act hemin40_mod o lect hemin40_mod o nO The detailed list of fortran routines The user will have to add the routines associated with his new grid The compilation instructions duplicated for each of th
10. horizontal high resolution a correction accounting for the differences in elevation between the climate model and the ice sheet model grids is necessary This correction can be performed off line as described in Section 5 or during run time In the case the correction is performed off line the user may set filel NH_40km_GRISLI xyz file2 NH_40km_GRISLI xyz Therefore no climate correction of the climate field will be performed during the initialisation since the correction will be equal to zero Following Equation 4 in which So would correspond to file2 and S t would correspond to file1 for t 0 A S t So 0 47 Therefore Teor t To 48 with T corresponding to the initial input downscaled surface temperature On the contrary if the correction is not performed off line the user will have to provide the climate model topography interpolated on the GRISLI grid domain and set Centro Euro Mediterraneo sui Cambiamenti Climatici CO Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers filel NH_40km_GCM on GRISLI xyz file2 NH_40km_GRISLI xyz Since the two topographies are different a climate correction will be performed during the initialisa tion S t So 0 49 Therefore Teor t To A S t So 50 Basal water amp eaubasalel nom du premier bloc eau basale coulement_eau T hwatermax 5000 000 infiltr 1 0000001E 03 ecoulem
11. is still being tested and might not work properly Reference climate forcing Centro Euro Mediterraneo sui Cambiamenti Climatici y Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers amp clim_act filel Snapshots GCM clim_B140_topo2_NH_40km xyz reference climate annual false type_precip tru filel reference climate annual true if climate file contains Tann Tjja type_precip true if liquid precip false snow The climate forcing i e surface temperature and precipitation is specified in this block The structure of the forcing file has been described in Section 5 annual if set to false the climate forcing file has to have monthly values 24 columns 12 for precipitation and 12 for temperatures In this case the annual mean is computed during run time If set to true the file should have mean annual values 1 column for precipitation and 1 column for temperature E type_precip corresponds to the nature of precipitation If set to true the prescribed precipitation corresponds to the total liquid precipitation while if set to false the routine assumes that the user provided the snowfall Transient climate forcing amp clim_forcage nom du bloc forcage_filel clim anoml 21k 0k dat mincoefbmelt 0 maxcoefbmelt 2 rapbmshelf 5 filforc b140_topol_150 dat tempgrad 0 005 tempgrjul 0 004 forcage_filel sn
12. resolution dx 40 km and the sea level here 0 m the third line is a header line where MKO correspond to the land sea mask and has to bet set to 0 everywhere for all Northern Hemisphere simulations rap precip corresponds to the precipitation ratio between the reference precipitation and the snapshot prescribed in this block deltaTann is the mean annual anomaly between the snapshot and the reference climate and deltaTjuly is the July temperature anomaly between the snapshot and the reference climate For example in the case of a transient simulation performed between the Last Glacial Maximum 21000 years and today the precipitation ratio will be Prem APLem gt 55 Ok the temperature anomalies will be computed as ATrom Trem Tor AlSzam Sok 56 AT Ey Trak Ti Stem Sor 57 Note that they both account for elevation difference between the two snapshot During run time the temperature is then modulated by the climate index T t Tok t A S t Sok ATiam X climindex t 58 Centro Euro Mediterraneo sui Cambiamenti Climatici CO Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers Similarly precipitation is interpolated and modulated by the same climate index The climate index must be computed such as for t 0 climindex 1 therefore T to Tram and for t n climindex 0 therefore T tn Tor A S t Sox The structure of the index file sh
13. restarting from T and Hwat 3 gt T only reprcptr name of the restart file This block is the first block of the whole namelist and includes all the parameters to be set related to the name of the simulation and the restart process E runname name of the simulation MUST HAVE 8 characters at maximum m icompteur switch to restart or not from a pre existing restart file reprcptr name of the restart file The restart files have the extension cptr and are written during run time in GRISLI Fichier CPTR at a frequency that have to set manually in the routines GRISLI SOURCES out_cptr_mod F90 To restart from one of those files the user may set the name of the restart file in this variable and change the cptr extension to CPTR For example if a previous simulation called NHtestO was run for 30 kyrs then the last restart file produced during run time will be NHtest0 k030 cptr To restart from this file the user may set icompteur 1 or 2 or 3 reprcptr NHtest0 k030 CPTR then the user may change the extension of the restart file to gt cd GRISLI Fichier CPTR gt mv NHtest0 k030 cptr NHtest0 k030 CPTR GRenoble Ice Shelf and Land Ice model a practical user guide In addition the user will have to set the initial model time in years to the time written in the first line of the restart file In the example above this initial model time will be set to 30000 See the explanations for the timesteps namelist block
14. the TEMPS NETCDF dat and TEMPS HZ dat located under GRISLI SOURCES Fichiers parametres fin des commentaires number of classes to be written 2 3 4 5 class see variables class in LISTE_VAR_NCDF dat 20 ndtsortie lecture de dtncdf de sortie 1000 2000 25 35 65 npredeft lecture de predef_tsort Centro Euro Mediterraneo sui Cambiamenti Climatici CO Centro Euro Mediterraneo sui Cambiamenti Climatici N CMCC Research Papers Three sections have to be parametrised 1 variable class the variables in the file LISTE VAR NETCDF dat have been classified according to their dimensions in five distinct classes 1 Var 2D to be output regularly 2 Var 2D to be output regularly and in case of debugging 3 Var 3D to be output regularly 4 Var of debugging only 5 Var output only once at the beginning and once at the end of the simulations In the example above 2 classes will be output i e class 1 and class 2 The order of the classes is not important and if the user wants to output class 3 and class 4 then this section will look like fin des commentaires 2 number of classes to be written 3 4 5 L 2 class see variables class in LISTE_VAR_NCDF dat ndtsortie Corresponds to the regular frequency of output given in years In the example above three frequencies have been set 20 years 1000 years and 2000 years Those frequencies can be assigned to the variables in th
15. the framework of the climate service in contract with the Swedish Nuclear Fuel and Waste and Management Company The following document provides a detailed description of the model and some indications to set up and run simulations Table of Contents NE gt e a ID AA E IR AAA A S 3 User GUE o a a i d a aiana D a E a A Aa a Mo o aa de dd de 4 1 Before installing asss ssamen eR OU EE RAPER eS 4 2 GRISLI directory structure 24 4 4 ee 5 3 Compiling and TUNNING s oi st a ak es 5 4 GRISLl grid demain s 244 ana daa a ee 8 5 GRISLI initial GONGHONE lt 4 4 caca aa a a e e eed de 9 6 GRISLI boundary conditions o oo a 13 7 GRISLI namelist example 13 8 Output variables and output frequency 18 Numeric and theory 44 21 Tutorial setting a Northern Hemisphere simulation 33 PHYSICS ODIOS 22 44 be us a ee 33 OMPI e a a eee Ra eee da te Dora ee 35 Namelist settings o 2 du da La due ou pu be die a 35 Tutorial how to create a new grid 53 o Gee E a a Pe eras no Sy a ere ee KORO ea ote eee e 56 o A i Ge wee eee ae od Da poesie 58 GRenoble Ice Shelf and Land Ice model a practical user guide Preamble GRISLI stand for GRenoble Ice Shelf and Land lce model It has been developed by Ritz et al 2001 at LGGE in Grenoble It is able
16. under GRISLI INPUT Forcage At each time step the depth of the base of the ice shelf is known Hraselt H i j t x Z 52 which calculates the ice thickness H i j t below the sea level based on Archimedes s law accounting for floatation with p as ice density and p as water density The temperature at the corresponding depth is interpolated between the ocean vertical levels Hbase i Js t ocdepth z ocdepth z 1 ocdepth z 99 Ti i j t Toe z Toelz 1 Toc z where ocdepth corresponds to the depth of each of the ocean levels prescribed in the input vertical temperature profile The melt rate is calculated as follow coyrlul Tur Ti t 4 Lex T t Tr oa m t where u takes an averaged value for the area of interest m s Tr is also set with an averaged value at a depth of interest yr is set to the ad hoc value of 1074 following literature values For example in the Arctic ocean there is no mixed layer as in the North Atlantic or around Antarctica therefore an averaged temperature below the halocline at 500m seems reasonable bmelt_type 3 Same as for bmelt_type 2 but using three dimensional vertical temperature maps one for each of the ocean level i e T lon lat z The file has the same structure as for the previous option but has n column containing temperature for each ocean vertical level The file is located under GRISLI INPUT Forcage This option
17. 0 tbegin 0 dtmin 2 e 2 dtmax Iss dtt 5 la testdiag 0 02 0 016 4 tous les temps en annees tbegin et tend debut et fin du run pour equation masse pas de temps mini gt dtmin maxi gt dtmax dtt pas de temps long testdiag pour g rer le pas de temps dynamique dt ordres de grandeur a moduler selon dx 40 km dtmin 2 e 2 dtmax 1 dtt 5 tesdiag 0 02 In this blocks the user determines the length of the run and set the time steps related to the resolution of the grid to satisfy the CFL conditions tend corresponds to the length or the run expressed in years If the user needs to perform a run of 100 kyrs then tend 100000 Alternatively the model also handles real time e for example if the user wants to perform a transient simulation based on input data evolving from the Last Glacial Maximum until today tend 0 and tbegin 21000 tbegin same as above For a 100 kyrs run long begin 0 m dtmin corresponds to the minimum time step needed by the mass equation and expressed in year At 40 km dtmin 2 e 2 dtmax maximum time step required by the mass equation At 40 km dtmax 1 i e one year dtt corresponds to the time step at which all the variables will be calculated and at which all the routines will be called in the main fortran routine main3D 0 4 40km F90 It can be different from dtmax For example if tbegin 0 the
18. 00000 1060000 0000000000 107 55830413905845 58 100928066994435 GRenoble Ice Shelf and Land Ice model a practical user guide The file contains no header and the red letters are only indicated here for a better understanding i and j correspond to the number of points in the X and Y directions respectively X and Y correspond to the coordinates of the grid expressed in native Cartesian distances km in the Lambert Equal Area projection Lon and Lat correspond to the longitudes and latitudes of the Cartesian coordinates X Y in the Lambert Equal Area projection In the vertical both the domain has 21 levels evenly spaced along the Z axis scaled to the ice thickness Furthermore there are 4 vertical levels in the bedrock layer 150 W 180 150 150 W 180 150 E 120 E ASAS APS SR Bedrock depression m Surface elevation m Ice thickness m Figure 2 Example of projected initial topography over the Eurasian domain at 20 km resolution bedrock topography surface elevation and ice thickness 5 GRISLI initial conditions The following initial conditions i e topography and climate forcing have to be projected on the horizontal grid before running GRISLI does not directly handle the projection procedure Initial topography GRISLI needs three initial topographic conditions namely surface topography bedrock topography i e topography below the existing ice sheets over ice free regions bedrock e
19. 241 84 9069672 where the first column corresponds to the i index coordinate number the second column corre sponds to the j index coordinates number and the third column holds the GHF values given in mW yr A second boundary condition used by GRISLI is a map of present day sediment thickness taken from Laske and Masters 1997 The original file is available at http igopweb ucsd edu gabi sediment html The data are provided on a grid of 1 x 1 hori zontal resolution and have to be interpolated on the GRISLI grid before running Figure 5a The sediment map is used by GRISLI to determine the possible areas where ce streams might occur The sediment map is related to the basal hydrology conditions see the section about Basal Hydrology below 7 GRISLI namelist example To each of the grid implemented in GRISLI corresponds a namelist encompassing all the parameters that the user may change An example of namelist is displayed below See the Tutorial section for the explanation of each namelist block Centro Euro Mediterraneo sui Cambiamenti Climatici OO Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers 150 W 180 150 E 120 E 105 E 120 W 90 W 90 E 60 W 75 E 30 W 60 E 45 E 9 Pe FF LPF POS SS OH LH SS FSF ES Sediment thickness m Figure 5 Map of sediment thickness from Laske and Masters 1997 interpolated on the Eurasian 20 km grid
20. CESM which is turned into snow using a density of 917 g cm In GRISLI ablation is parametrised using the Positive Degree Day PDD semi empirical method Reeh 1991 This method is based on an empirical relationship between the number of positive degree days computed from annual mean and July surface air temperature and the snow and ice melting rates which depend on the melting factors Csnow and Cice derived from observations The number of positive degree days PDD is given by 1 gt T Ty PDD I es EA dTdt 45 Ooy 2T lyr JO 20 where Ty is the daily temperature and o the standard deviation of the daily temperature usually set to 5 C This formulation allows for positive temperatures even when the average daily temperature is below the melting point As a consequence it may also overestimate melting in general Bougamont et al 2007 The daily temperature T4 is reconstructed from annual mean and July temperatures Tann and Truy by assuming that the annual temperature cycle follows a cosine function Figure 7 Talt T Ann TJuly Tann cos 2 7 t 365 46 Ablation is then calculated in several steps GRenoble Ice Shelf and Land Ice model a practical user guide TEMPERATURE Figure 7 Reconstructed annual temperature cycle used to calculate the number of positive degree days The annual cycle is retrieved from the annual mean temperature and the July temperature Figure from Reeh 1991
21. CMCC Centro Euro Mediterraneo sui Cambiamenti Climatici Research Papers Issue RP0249 January 2015 Ocean Modeling and Data assimilation Division By Florence Colleoni Centro Euro Mediterraneo sui i Cambiamenti Climatici Bologna Italy flocolleoni gmail com We thank Catherine Ritz for providing the GRISLI code and for the helpful discussions and Claudia Wekerle on the description of the model physics We gratefully acknowledge the support of Italian Ministry of Education University and Research and Ministry for Environment Land and Sea through the project GEMINA GRenoble Ice Shelf and Land lce model a practical user guide SUMMARY GRISLI has been developed by Ritz et al 2001 at LGGE CNRS France It is able to simulate both grounded and floating ice The grounded part uses the Shallow Ice Approximation SIA Hutter 1983 whereas ice shelves and ice streams are simulated following the Shallow Shelf Approximation SSA MacAyeal 1989 The ice shelf formulation in GRISLI allows for a more realistic calculation of the ice sheet growth and particularly of the advance of ice onto the shallow continental shelves in both hemispheres high latitudes GRISLI has been validated over Antarctica Greenland and successfully applied to study the inception of ice sheets during the last glacial cycle At CMCC GRISLI has been used to simulate the glacial inceptions of the last two glacial cycles and the penultimate glaciation in
22. EAS Ice streams and ice shelves are characterised by fast flow and low surface slopes The sliding velocity eq 22 depends strongly on the surface slope and thus fast flowing areas would not be represented well in the model For that reason we treat ice stream and ice shelf areas separately in the model by applying the shallow shelf approximation SSA MacAyeal 1989 Those regions are called fast flowing areas Due to the fact that we use a horizontal resolution of 20 km and 40 km in our model setup GRISLI cannot resolve individual ice streams Thus a grid point in GRISLI is considered as ice stream grid point if it has the large scale characteristics of an ice stream area With the SSA the horizontal velocity is independent from depth which is reasonable for ice shelves and fast flowing ice stream areas This leads to the following elliptic equation for horizontal vertically averaged velocity u Ou Ov _ Ou Ov Os 3x 27n e By fan e poli Th s 24 2 v u a dv du as ae 22 at A a8 2 Oy 7 Oy a i Ox Ox pg Oy Tby 25 where 7 is the ice viscosity vertically averaged over the ice thickness 7 and 7 stand for the basal drag This set of equations is used both for ice stream areas and ice shelves Ice shelves In case of ice shelves the basal drag m is set to zero since ice shelves are floating on the water without any resistance In every time step the position of i
23. F_INC and NCDF_LIB respectively such that NCDF_INC usr local netcdf 4 1 3_ifort include NCDF_LIB usr local netcdf 4 1 3_ifort lib STEP 2 Compiler setting the name of the compiler in the variable IFORT IFORT gfortran for example Centro Euro Mediterraneo sui Cambiamenti Climatici O1 Centro Euro Mediterraneo sui Cambiamenti Climatici D CMCC Research Papers GRISLI SOURCES BLAS Eurasie20_files Fichiers parametres Hemin40_files Netcdf routines New remplimat INPUT HEMIN 40 EURASIE20 Forcag e RESULTATS L Netcdf Resu Param bin Figure 1 Description of GRISLI main directories STEP 3 Compilation options setting the options to compile the model according to your com piler and your operative system in the four variables ARITHM FT LK F_NETCDF For example to compile GRISLI using the IBM Fortran and C compilers on 64 bits Macintosh architecture ARITHM EJE LK F_NI ETCDF 02 xHost IFORI l l Sl IFORI IFORI l ARIT ARIT ARIT no prec div mk1 m64 THM c IS NCDF_INC q THM g THM c IS NCDF_INC Note that in the ARITHM variable we include mk1 option corresponding to the famous IBM mathematical MKL libraries provided with the IBM Fortran compiler Those libraries contains the BLAS routines which are necessary to run GRISLI To compile GRISLI using the GNU Fortran and C
24. J mol Pa Pa Pa yr Pa Pa m s m yr Pa 910 kg m73 1028 kg m73 1000 kg m73 9 81 ms 8 3145 J mol C 1 Centro Euro Mediterraneo sui Cambiamenti Climatici SJ Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers References Alvarez Solas J Montoya M Ritz C Ramstein G Charbit S Dumas C Nisancioglu K Dokken T and Ganopolski A 2011 Heinrich event 1 an example of dynamical ice sheet reaction to oceanic changes Clim Past 7 1297 1306 Bougamont M Bamber J L Ridley J K Gladstone R M Greuell W Hanna E Payne A J and Rutt 2007 Impact of model physics on estimating the surface mass balance of the Greenland ice sheet Geophys Res Lett 34 L17501 Charbit S Ritz C and Ramstein G 2002 Simulations of North ern Hemisphere ice sheet retreat sensitivity to physical mecha nisms involved during the Last Deglaciation Quaternary Science Reviews 21 243 265 Dumas C 2002 Mod lisation de l volution de l Antarctique depuis le dernier cycle glaciaire interglaciaire jusqu au futur importance relative des diff rents processus physiques et r le des donn es d entr e PhD thesis Universit Joseph Fourier Grenoble 1 Greve R and Blatter H 2009 Dynamics of Ice Sheets and Glaciers Springer Verlag Berlin Holland P R Jenkins A and Holland D M 2008 The response of ice shelf basal melting to var
25. al resolution in km Climate forcing The ice sheet model is forced by steady state monthly air temperatures and monthly total precipi tation water equivalent from any climate data provider observations simulations Both climate fields have to be provided for each of the 12 month of the year Climatic means when needed by GRISLI are performed during run time Temperatures have to be provided in degree Celsius Precipitations have to be provided in meter year in water equivalent The conversion into snow and ice is performed during run time using the appropriate density Both climate fields have to be projected and interpolated on the GRISLI horizontal rectangular grid before running Both variables have to be written in an ASCII files containing 3 headers and 24 columns 12 column for precipitations first followed by 12 columns for temperatures GRISLI reads the file following the Fortran format 24 e12 5 1x 140000 000 Hemisphere Nord climate forcing CESM fv09 B_case 56700 24 210 270 20 000000000000000 Monthly precip m an and monthly T2M C 2 066 1 887 he 2 181 5 430 4 537 ee 6 504 GRenoble Ice Shelf and Land Ice model a practical user guide The climate forcing generally comes on a standard Mercator longitude latitude grid at a resolution varying from 0 5 to 1 The original climate forcing grid resolution is generally coarser than the horizontal resolution of the ice sheet model As a consequence the
26. amp eaubasalel nom du premier bloc eau basale ecoulement_eau T hwatermax 5000 000 infiltr 1 0000001E 03 coulement eau false gt modele bucket sinon equ diffusion hwatermax hauteur d eau basale maximum dans le sediment m infiltr est la quantite d eau qui peut s infiltrer dans le sol 1 amp param_hydr nom du bloc parametres hydrauliques hmax_till 20 00000 poro_till 0 5000000 kondo 1 000000E 06 hmax_till m epaisseur max du sediment poro_till porosite du sediment conductivite du sediment kond0 m s amp drag_hwat_cont nom du bloc dragging hwater contigu hwatstream 50 cf 1 e 5 1 e 4 tobmax 1000 toblim 0 7e5 0 25e5 hwatstream m critere de passage en stream en partant de la cote si hwater tobmax gt hwatstream Pa frottement maxi f cf coefficient de la loi de frottement fonction Neff 1 1 toblim Pa tes pour les iles amp drag_hwat_s edim nom du bloc dragging hwater sediment filesedim nosediment_ij_hemin40 dat hwatstream 50 cf 2 e 5 1 e 4 tobmax 1000 toblim 0 7e5 0 25e5 seuil_hwater 250 a practical user guide ici Centro Euro Mediterraneo sui Cambiamenti Climat ici Centro Euro Mediterraneo sui Cambiamenti Climat CMCC Research Papers seuil_sedim 150 seuil_neff 350 e5 hwatstream m critere de passage en stream en partant de la cote si hwater gt hw
27. apshot climat 1 mincoefbmelt borne mini de coefbmshelf maxcoefbmelt borne maxi de coefbmshelf filforc fichier de forcage tempgrad lapse rate annuel tempgrjul lapse rate d ete GRenoble Ice Shelf and Land Ice model a practical user guide In GRISLI it is possible to perform steady state and transient simulations The choice of clim_forcage in the module_choix leads to transient simulations Transient simulations consist in interpolating the climate forcing between two snapshots modulated by a climate index varying in time 4190 for example or sea level etc In this case the user may specify one or several climate snapshots representative of two different periods between which to interpolate the climate To use this module the forcing file might be computed as an anomaly relative to the prescribed reference climate i e the climate forcing file set in the clim_act block Therefore the structure of the forcing file should be 21000 00 H Nord forcage climatique k115 k125 CESM 58081 4 241 241 40 0 0 MKO rap precip deltaTann deltaTjuly 0 0 487772 0 2986 0 2491 0 0 467371 0 2968 0 2527 where the first line should states the date of the snapshot here 21 000 years real paleotime and the rest of the line is a comment about the climate forcing the second line should states the total number of lines in the file NX NY the number of column 4 the grid domain dimensions e NX and NY the horizontal
28. areas the grounded ice is generally thinner than Hcoup AS a con sequence when the ice starts to float some calving occurs As the horizontal velocities increase substantially in the floating part the ice flux is too large which does not allow the ice to reach a thickness larger than H p Therefore an additional test is performed in that case and the testing procedure is applied only if none of the neighbour point is grounded GRenoble Ice Shelf and Land Ice model a practical user guide At last if some nodes in the middle of the ice shelves have an ice thickness lower than H p Some polynias develop locally Those polynias can grow and divide the ice shelf This issue is avoided by applying calving only if a point is located at the front of the ice shelf In GRISLI there are 5 different options to treat the calving For example one of the option is to test whether or not the ice thickness in the upstream point is larger than Hou or the mass balance bm bmelt is positive Those options are described more into details in the Tutorial section BASAL MELTING UNDER THE SHELVES Because they are floating over oceans the ice shelves are naturally affected by the ocean warmth at their base The basal melting is generally important near the grounded line area because the depth of the ice shelves is large and therefore reaches area were ocean temperature is more elevated than at the surface The fresh water resulting from the ice melting i
29. ariable L dd t isortie dtncdf interv l ice thickness GRenoble Ice Shelf and Land Ice model a practical user guide where the frequency of output is determined according the three flags isortie dtncdf interv isortie determines whether or not the variable will be output If isortie 1 gt output if isortie 0 gt no output dtncdf indicates the periodicity of output for the variables For example dtncdf 1 means that the variables will be output at the first time step frequency ndtsortie indicated in the TEMPS NETCDF dat dtncdf 2 means that the variables will be output at the second time step frequency ndtsortie indicated in the TEMPS NETCDF dat dtncdf 0 means that the output will occur only on the intermediate time steps npredeft indicated in the TEMPS NETCDF dat interv indicates the time interval at which outputs are written interv 1 the outputs are written at the first and last time steps only interv 0 the outputs are written at the first time step only interv 1 the outputs are written at all time steps specified in the file TEMPS NETCDF dat For example 1 0 outputs written for the first time step and all intermediate time steps specified in TEMPS NETCDF dat 1 outputs written at the first and last time step 1 no output written 0 outputs written at the first time step only 1 outputs written at all time steps first intermediate and last RO The time step frequency is defined in
30. atstream uf 1 cf coefficient de la loi de frottement fonction Neff tobmax Pa frottement maxi 1 1 toblim Pa tes pour les iles seuil_hwater m seuil hwater pour avoir glissement sur zone sediment seuil_sedim m seuil epaisseur sediment pour avoir glissement sur zone sediment amp drag_jorge nom du bloc dragging_jorge filesedim nosediment_ij_hemin40 dat hwatstream DU cf 2 e 5 1 e 4 tobmax 1000 toblim 0 7e5 0 7e5 seuil_sedim Ty 0 0001 hwatstream m critere de passage en stream en partant de la cote si hwater gt hwatstream cf coefficient de la loi de frottement fonction Neff 1 1 tobmax Pa frottement maxi toblim Pa tes pour les iles amp ghflux geothermal heat fluxes fileghf ijphi_hemin40 dat coefghf 1 fileghf input geothermal heat fluxes map coefghf to modulate the amplitude of the geothermal heat fluxes 1 amp calving l nom du bloc calving m thode Vincent Hcoup 200 ifrange 3 meth_Hcoup 0 Hcoup epaisseur de coupure ifrange 0 gt NO CALVING ifrange 1 gt traitement de Vincent avec ice shelves frangeants biased ifrange 2 gt ice shelves frangeant seulement si bm bmelt positif biased ifrange 3 gt option 2 corrected but no test on ice thickness of upstream point ifrange 4 gt option 3 but with test on upstream point ice thickness point meth_Hcoup gt Hcoup depends on
31. ature fields have to be corrected for these elevation changes by means of a lapse rate value downscaling TisM cor Tecm on 1smM A Seom on 1sM Sism Q where Tism cor and Toom on rsm are the final downscaled temperature field and the inter polated temperature field respectively S7sy1 denotes the elevation of topography on the high Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers resolution ice sheet model grid and Sozsm on 1sm is the elevation of the CESM topogra phy interpolated on the high resolution ice sheet model grid At last A corresponds to the atmospheric lapse rate C km used to correct temperature for elevation difference The same procedure can be used to correct precipitations The correction can be performed using various methods Here we apply the same method as GRISLI uses during run time based on Charbit et al 2002 PisM cor Pacm on 1sm exp y Tecm on 1sm TISM cor 3 where Por and Po stand for the corrected and initial precipitation respectively y corresponds to the precipitation correction factor Note that the use of exponential function in the precipitation formulation eq 5 is motivated by the saturation pressure of water vapour in the atmosphere Clausius Clapeyron relationship which increases roughly exponentially with temperature 150 W 180 150 E 120 E 105 E Centro Euro Mediterraneo sui Cambiamenti Climatici 45 E w
32. bove the pressure melting point and the effective pressure is high we choose a constant value of Ky 1078 m s If the effective pressure is low the hydraulic conductivity increases K if N gt 100b 0 i gt 100 bar 31 Ko 1002 if N lt 100 bar The variation of the hydraulic head with time is governed by the horizontal variation of the water flux and by two source terms the basal melting rate bme computed from the heat equation eq 8 and the infiltration rate into the soil Tin piter Ohy OVex oV a alas a F Tin pittr 32 ICE TEMPERATURE The ice temperature is required to calculate the flow law coefficient Bar the basal melting rate F below the grounded regions and the basal ice velocity us since sliding occurs only in areas where the basal ice temperature is at the melting point The time dependent heat equation 8 is solved in the ice sheet and accounts for horizontal and vertical advection of ice vertical diffusion and production of heat through deformation At the surface the ice temperature is determined by climatic conditions which also depend on surface elevation The ice temperature 10 m below the ice surface is almost similar to the annual mean surface air temperature Tann Thus we use Tann as a surface boundary condition for ice temperature At the base of the ice sheet a geothermal heat flux GHF is prescribed At the ice bedrock interface we distinguish between warm and cold conditions B
33. by the user The projection informations are not needed by GRISLI but are useful in case the user wants to reproject the outputs from the simulations on a standard longitude latitude map The number of points in X and Y directions depends on the size of the domain and on the horizontal resolution Those numbers are decided by the user and have to be provided to the model see Tutorial section The X Y coordinates of the grid must be provided in a separate coordinate ASCII file in native Cartesian distances km also containing the index number of each pixel i j as well as the corresponding projected longitude and latitude Lon Lat For example the GRISLI grid covering part of Northern Eurasia projected into Lambert Equal Area geographical projection centre of projection 0 90 N Figure 2 for an example at 20 km resolution contains 210 in the X direction and 270 points in the Y direction The corresponding coordinate file located under GRISLI INPUT EURASIE20 has to be structured as follows i j X Y Lon Lat 1 1 830000 00000000000 4320000 0000000000 349 12431166414439 49 734095498049626 2 1 810000 00000000000 4320000 0000000000 349 38034472384487 49 769806622114523 3 1 790000 00000000000 4320000 0000000000 349 63680754375133 49 804671832343111 4 l 770000 00000000000 4320000 0000000000 349 89369090450970 49 838689334153315 5 1 750000 00000000000 4320000 0000000000 350 15098548888295 49 871857370085763 210 270 3350000 00000
34. ce shelves is determined with a flotation criterion which is based on Archimedes principle of floating bodies pw SL B pH 26 where p and pu are the ice and water density respectively SZ and B stand for sea level and the base of the ice sheet respectively Calving of the ice shelf is determined by a thickness criterion As soon as the ice thickness in an ice shelf node at the front reaches a certain threshold H ai the node is cut off We set Heaw to 200 m which corresponds to values observed at the large ice shelves of Antarctica The oceanic basal meting rate below the ice shelf bme is determined by the temperature of the ocean In this version of GRISLI we use a simplified approach assuming a melting rate which only depends on the water depth Above 450 m bu is set to 0 2 m yr while below 450 m depth bme is set to 2 45 m yr Ice streams In case of ice streams the basal drag in equations 24 and 25 is described by a friction law It depends on the vertically averaged velocity u and the effective pressure N eq 23 scaled with the basal drag coefficient cy Th Cf NU 27 The position of the ice stream grid points is determined with the following three criteria from which at least one has to be fulfilled Centro Euro Mediterraneo sui Cambiamenti Climatici dl Centro Euro Mediterraneo sui Cambiamenti Climatici N CMCC Research Papers The grid point is located in a narrow valley Firs
35. compared to the vertical normal stress 7 This reduces the vertical momentum equation to a balance between the vertical gradient of r and the gravity force hydrostatic approximation Integration gives the hydrostatic pressure equation p pg H 2 14 In GRISLI like in other large scale ice sheet models the shallow ice approximation SIA Hutter 1983 is used The SIA is a simplification which is reasonable when the horizontal length scale is much larger than the ice thickness We assume that the flow regime is mainly characterised by bed parallel shear Thus the relevant components of the deviatoric stress tensor are the shear Centro Euro Mediterraneo sui Cambiamenti Climatici CD Centro Euro Mediterraneo sui Cambiamenti Climatici N N CMCC Research Papers stresses in the horizontal plane 7 and 7 Then the horizontal components of the momentum balance simplify to Tsz Op 0H dz r I Og gt OTyz _ Op OH Oz Oy e Oy 16 The deviatoric stress Tsy is related to the strain rate by a nonlinear Glen type flow law an equivalent expression holds for 7 and z 1 04 Ow ne Ery 2 5 a Bar Tpmp T cer 17 where Bar is the temperature dependent flow law coefficient n is the flow law exponent here we use n 3 7 is the effective shear stress given by Ae Tiz 18 The ice temperature Tm is the absolute temperature corrected for the dependence on the melting point on pressur
36. compilers on 64 bits Macintosh architecture GRenoble Ice Shelf and Land Ice model a practical user guide ARITHM O2 ffr line length non mtune nativ mfpmath sse FT IFORT ARITHM c I NCDF_INC LK IFORT S ARITHM F_NETCDF IFORT S ARITHM c I NCDF_INC Note that the MKL option is not specified In this case the user will have to follow STEP 4 STEP 4 BLAS library include or not the BLAS library in the compilation options If the user does not use the IBM Fortran and C compilers some of the mathematical functions needed to solve the elliptic matrix included within the MKL libraries will be missing Therefore the user will have to compile the GRISLI SOURCES BLAS routines provided with the code if the directory is missing please report to the previous sections The routines needed by GRISLI are already listed within the Makefile To compile them the user will have to manually remove all the x o files and execute the compilation At last the user will have to manually include the Liste BLAS libraries in the compilation lines of the selected grid do main of the Makefile For example to compile the routine to run over the Northern Hemisphere Hemin 40 Dim_hemin40 mod_communs mod_ell S Liste_hemin40 diagnoshelf Liste_Netcdf routines_communes S routine_elliptiques Liste_BLAS Centro Euro Mediterraneo sui Cambiamenti Climatici SJ S LK o RUND Hemin 40
37. depth bmelt_type 1 bmelt1 0 bmelt2 0 depthlimit 2000 ocn_profile TEMP_vert_Chukchi dat ocn_map Tocn_B140_topol_artic_20km_noice_GRISLI xyz bmelt_type 1 gt prescribed basal melting depending on depthlimit bmelt_type 2 gt basal melting calculated using an ocean temperature vertical profile bmelt_type 3 gt as for 2 but using bidimensional ocean vertical temperature bmelt1 basal melting above depthlimit bmelt2 basal melting below depthlimit depthlimit depth threshold to precribe basal melting ocn_profile ocean temperature vertical profile ocn_map vertical ocean temperature map 1 1 1 1 pl 1 1 1 amp bmelt_nor_reg1 nom du bloc lbmelt_Ross bmgrz_Ross lbmelt_FRis bmgrz_FRis ln lo 0 01 N amp clim_forcage nom du bloc forcage_filel forcage_file2 mincoefbmelt maxcoefbmelt rapbmshelf 5 filforc b140_topo1_150 dat tempgrad 0 0 tempgrjul 0 0 forcage_filel snapshot climat 1 mincoefbmelt borne mini de coefbmshelf maxcoefbmelt borne maxi de coefbmshelf filforc fichier de forcage tempgrad lapse rate annuel tempgrjul lapse rate d ete Il 0 2 ll amp clim_act filel Snapshots GCM clim_B140_topo2_NH_40km_withice_GRISLI_lapse5 xyz reference climate ci Centro Euro Mediterraneo sui Cambiamenti Climat Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Re
38. e Tomp T 9 76 107 pg S 2 19 The flow law coefficient B 47 follows an Arrhenius relationship Ea 1 1 Bar Baroexp z FJ 20 pmp where Baro is a coefficient E is the activation energy and R the universal gas constant Combining the above equations and integrating over z gives an equation for the horizontal velocity u 5 u p9 VS V8 VS Ey Bar Tomp S 2 de up 21 B E is a parameter the flow enhancement factor set to 3 in this model version We assume that sliding of ice over the bedrock occurs only if the basal ice temperature reaches the melting point The basal velocity ug is expressed as a function of basal shear stress and effective pressure N following Weertman 1957 up k pgH VS VS 2 VS N 1 22 where k is an adjustable parameter set to 5 1071 in this version of GRISLI The effective pressure N is given by N pgH pw 23 where p is the sub glacial water pressure Again equation 22 is only valid for regions where the basal ice temperature reaches the melting point If the basal ice temperature is below the melting point the basal velocity is set to zero In GRISLI the sliding law described above is only valid for small velocities lt 50 m yr Regions with larger velocities are considered as fast flowing areas and treated separately see next section GRenoble Ice Shelf and Land Ice model a practical user guide EQUATION OF MOTION FAST FLOWING AR
39. e 5 1 e 4 tobmax 1000 toblim 0 7e5 0 7e5 seuil_sedim Ls 0 0001 hwatstream m critere de passage en stream en partant de la cote si hwater gt hwatstream cf coefficient de la loi de frottement fonction Neff tobmax Pa frottement maxi toblim Pa tes pour les iles Those three blocks are related to the sediment and basal water thresholds needed to trigger basal dragging and in particular the ice streams i e the criteria to treat grounded ice with the SSA Each of them correspond to fortran routine located under GRISLI SOURCES Hemin40 files The choice of the dragging routine has to be made in the routine GRISLI SOURCES Hemin40 files module choix hemin40 0 4 f90 BEFORE compiling the model The first block drag_hwat_cont contains parameters to trigger ice streams based on the melt water level at the base of the ice sheet The two following blocks drag_hwat_sedim and drag_jorge account for sediment thickness at the base of the ice sheet and the amount of water contained in those sediments to trigger ice streams E filesedim is the sediment map interpolated on GRISLI grid and located under GRISLI INPUT HEMIN40 E hwatstream is the minimum water level above which ice streams are trigger if hwater gt hwatstream m cf if the dragging coefficient function of the effective pressure Nef f a tobmax maximum basal dragging in Pascal toblim maximum pressure for islands in Pascal E seuil_hwater
40. e density is constant and accounting for ablation and accumulation Surface Mass Balance SMB at the ice surface as well as melting at GRenoble Ice Shelf and Land Ice model a practical user guide the ice base F leads to OH Hu dt Ox O H v Oy SMB F 10 The fundamental contributions that govern the mass balance of the ice sheet are m The surface mass balance SMB accumulation minus ablation which depends on climatic conditions surface air temperature and precipitation The basal melting rate F which depends on the pressure melting point below the ice sheet and the geothermal heat flux Below ice shelves basal melting results from oceanic heat flux and is denoted by bc m Calving of the ice shelf which depends on a thickness criterion EQUATION OF MOTION GROUNDED ICE Velocities within the ice sheet are rather small and thus acceleration can be ignored Assuming typical values for an ice sheet geometry and typical flow velocities the ratio between acceleration and horizontal pressure gradient is 10 15 see Greve and Blatter 2009 Equation 7 can then be written as OTxx j OTxy A Tzz Ox Oy Oz GU OTya OTyy OTyz _ r y 0z 2 O z o z o 22 DE y Ta E pg 13 Ox Oy Oz We split up the stress tensor into a deviatoric part 7 and an isotropic pressure 1 Tij Tij g ze Tyy Tz Oi Tij TP We assume that the shear stresses 7 and r are small
41. e file LISTE VAR NCDF dat 111 the variable will be output each 20 years Another example 1 2 1 the variable will be output each 1000 years ndtpredeft corresponds to intermediate time steps given in years In the example above outputs will also occur at model year 25 model year 35 model year 65 and model year 85 GRenoble Ice Shelf and Land Ice model a practical user guide Numeric and theory GRISLI is able to simulate both grounded and floating ice The grounded part uses the Shallow Ice Approximation SIA Hutter 1983 whereas ice shelves and ice streams are simulated following the Shallow Shelf Approximation SSA MacAyeal 1989 The ice shelf formulation in GRISLI allows for a more realistic calculation of the ice sheet growth and particularly of the advance of ice onto the shallow continental shelves in both hemispheres high latitudes GRISLI has been validated over Antarctica Ritz et al 2001 and successfully applied to study the inception of ice sheets during the Early Weichselian period Peyaud et al 2007 Alvarez Solas et al 2011 In the following we give a detailed description of the model Notation of the main variables is listed in Table 1 and Table 2 DEFINITION OF ICE FLOW IN GRISLI In GRISLI three different regions of ice flow are considered ice sheets ice stream zones and ice shelves Ice sheets and ice stream zones are regions with grounded ice whereas ice shelves are regions with f
42. e grid domains implemented in the model The user has to create the same block for his own new grid Hemin 40 Dim_hemin40 mod_communs mod_ell Liste_hemin40 diagnoshelf Liste_Netcdf routines_communes routine_elliptiques Liste_BLAS S LK o RUND Hemin 40 Dim_hemin40 mod_communs mod_ell Liste_hemin40 diagnoshelf Liste_Netcdf routines_communes routine_elliptiques L NCDF_LIB lnetcdf L NCDF_LIB lnetcdff Liste_BLAS Centro Euro Mediterraneo sui Cambiamenti Climatici Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers Appendix Geometric and physical variables S x y t So x y SL B x y t Bo x y t H zx y t Ho x y t P x y 2 Tay 2 6 To x y 2 t Tpmp 2 y 2 t u x y z t u x y t up z y t Mass balance and Climate Abl x y t Acc x y t Cice Csnow P T Ut PDD SMB x y t Tann Tyul Ta bmelt cst AAnn AJuly Y oO Table 1 Notation used in this report Surface elevation Initial surface elevation Sea level Bedrock elevation Initial bedrock elevation Ice thickness H S B Initial ice thickness Hydrostatic pressure Ice temperature Basal ice temperature Ice temperature relative to the melting point of pressure Ice velocity Vertically averaged horizontal ice velocity Sliding velocity at the base of the ice Ablation rate Accumulation rate Melting factor of ice Melti
43. e located under GRISLI SOURCES Hemin40 files and GRISLI SOURCES Eurasie20 files Those di rectories contains the main information related to each grid in the paradim hemin40_mod f90 and paradim euras20 mod f90 module module _ nord implicit none character len 7 parameter geoplace hemin20 geoplace character len 200 parameter dirnameinp Users GRISLI INPUT HEMIN40 character len 200 parameter dirforcage Users GRISLI INPUT Forcage dimensionnement grilles integer parameter nx 241 dimension selon x indice 1 integer parameter ny 241 dimension selon y indice j integer parameter nlmin 141 min nx ny pour matrice elliptique integer parameter nz 21 dimension selon z indice k integer parameter nzm 4 nb de points verticaux dans le socle integer parameter nn 241 pour des sorties lcdc integer parameter nxl nx 1 nyl ny 1 nzl nz 1 Centro Euro Mediterraneo sui Cambiamenti Climatici GO Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers pas horizontal 1 real parameter dx 40000 0 real parameter dy 40000 0 pas en X et Y enm end module module_nord The user has to adapt the various dimensions name and repository to its new grid Other routines such as the ones reading the input grid input geothermal heat fluxes file the sediment map and the climate forcing has to be edited The user may not forget to also
44. e sheet over the last 420 000 years Im plications for altitude changes in the Vostok region J Geophy Res 106 D23 Shapiro N and Ritzwoller M 2004 Inferring surface heat flux distributions guided by a global seismic model particular appli cation to Antarctica Earth and Planetary Science Letters 223 213 224 Turcotte D L and Schubert G 2002 Geodynamics Cambridge Univ Press New York Weertman J 1957 On the sliding of glaciers CJ Glaciol 3 33 38 GRenoble Ice Shelf and Land Ice model a practical user guide Centro Euro Mediterraneo sui Cambiamenti Climatici 2015 Visit www cmcc it for information on our activities and publications The Euro Mediteranean Centre on Climate Change is a Ltd Company with its registered office and administration in Lecce and local units in Bologna Venice Capua Sassari Viterbo Benevento and Milan The society doesn t pursue profitable ends and aims to realize and manage the Centre its promotion and Centro Euro Mediterraneo research coordination and different scientific and applied activities in the field of climate change study sul Cambiamenti Climatici Centro Euro Mediterraneo sui Cambiamenti Climatici CO
45. e user might check the installation of the following tools A fortran compiler GRISLI has been developed originally using the IBM fortran and C compilers namely ifort and icc GRISLI has been successfully tested with the open source GNU Fortran and C compiler gfortran and gcc The results produced by both compilers are almost identical when using the appropriate compiling and optimisation options Basic Linear Algebra Subprograms BLAS fortran libraries those libraries are needed to perform some of the mathematical calculations needed by GRISLI While this libraries came originally along with the IBM compilers an open source version available for download here http www netlib org blas The routines have been developed in Fortran 77 The open source libraries might not be included within GRISLI code and the user might have to download it NetCDF fortran libraries As GRISLI outputs are produced both in ASCII and NetCDF format the user should compile the NetCDF Fortran and C libraries available for download on the main NetCDF website with their own Fortran compiler Note that the NetCDF releases after 4 2 does not include Fortran and C libraries in the same package which have to be installed separately GRISLI runs under UNIX environment and has been successfully run on Macintosh operative systems GRISLI runs in serial on one single processor and has not been optimised to use multi threads processes To post process a
46. elow the melting point the heat flux at the interface consists of geothermal heat which is transmitted by vertical diffusion into the ice At the melting point the melting rate depends on the difference between the heat fluxes on both sides including heat generated by friction during the sliding If basal temperatures T are at the melting point they are held constant T Tomp if Tp gt Tomp 33 and the excess in heat is used for basal melting F This basal melting rate F is one of the components of the mass conservation equation eq 6 CALVING The spatial extent of ice shelves is controlled by the calving process In GRISLI the position of the ice shelf front is based on a thickness criterion Heoup The ice shelf front corresponds to the pixels where the ice thickness at the ice ocean boundary equals 0 A grid node is considered as part of the front if all the three following conditions are fulfilled floating Centro Euro Mediterraneo sui Cambiamenti Climatici y Centro Euro Mediterraneo sui Cambiamenti Climatici N CMCC Research Papers the ice thickness is not equal to zero if a neighbour point has an ice thickness equals to zero and is located in the ocean domain The method implemented by Peyaud 2006 consists in determining if a new ice shelf point should be cut or should thicken and become part of the ice shelf H gt H oup Ih this node is not calved it will be considered as the new ice shelf fron
47. ent eau false gt modele bucket sinon equ diffusion hwatermax hauteur d eau basale maximum dans le sediment m infiltr est la quantite d eau qui peut s infiltrer dans le sol m an 1 Those parameters control the amount of basal water at the base of the ice sheet ecoulement_eau is a logical If set to false the bucket model is used otherwise the water flow is handle by diffusion E hwatermax corresponds the maximum water level contained in the sediments sediment map input also in the namelist in meters infiltr is the maximum amount of water able to penetrate in the ground in m yr Hydraulic parameters amp param_hydr nom du bloc parametres hydrauliques hmax_till 20 00000 poro_till 0 5000000 kondo 1 000000E 06 GRenoble Ice Shelf and Land Ice model a practical user guide hmax_till m max sediment thickness poro_till sediment porosity kond0 sediment conductivity m s 1 This block set the hydraulic parameters used to drain sub glacial waters The basal hydrology in GRISLI is based on a Darcy type flow law stating that water flows through a porous medium till sediments proportional to a hydraulic pressure gradient This is why the properties of this till sediment layer at the base of the ice sheet have to be set m hmax_ti11 is the maximum sediment thickness considered in the drainage of the sub glacial waters m poro_till is the porosity of this sedim
48. ent layer at the base This can be varied if the user changes the nature of the till kondo is the sediment conductivity depending on the nature of the till Basal dragging and ice streams amp drag_hwat_cont nom du bloc dragging hwater contigu hwatstream 50 cf 1 e 5 1 e 4 tobmax 1000 toblim 0 7e5 0 25e5 hwatstream m critere de passage en stream en partant de la cote si hwater gt hwatstream cf coefficient de la loi de frottement fonction Neff tobmax Pa frottement maxi toblim Pa tes pour les iles amp drag_hwat_sedim nom du bloc dragging hwater sediment filesedim nosediment_ij_hemin40 dat hwatstream 50 GE 2 e 5 1 e 4 tobmax 1000 toblim 0 7e5 0 25e5 seuil_hwater 250 seuil_sedim 150 seuil _neff 350 e5 hwatstream m critere de passage en stream en partant de la cote Centro Euro Mediterraneo sui Cambiamenti Climatici Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers si hwater gt hwatstream cf coefficient de la loi de frottement fonction Neff tobmax Pa frottement maxi toblim Pa tes pour les iles seuil_hwater m seuil hwater pour avoir glissement sur zone sediment seuil_sedim m seuil epaisseur sediment pour avoir glissement sur zone sediment amp drag_jorge nom du bloc dragging_jorge filesedim nosediment_ij_hemin40 dat hwatstream 50s cf 2
49. et Q cold Q1 et Q warm 02 loi de deformation 2 amp loidef_2 exposant_2 1 temp_trans_2 10 enhanc_fact_2 10 coef_cold_2 8 313 Q_cold_2 4 0001 coef_warm_2 8 313 Q_warm_2 6 0001 A AH A F E Se OO HS 00 Hs 00 module deformation_mod_2lois exponent glen transition temperature ttrans enhancement factor sf GRenoble Ice Shelf and Land Ice model a practical user guide for temperatures lower than Temp_trans coef_cold Bat1 et Q cold 01 for temperatures greater than Temp_trans coef_ warm Bat2 et Q warm 02 The deformation block handles the parameters of the Glen flow law see Numeric and theory section applied to the grounded ice SIA In this version of GRISLI the flow law uses two sets of parameters reported in the namelist The use of the first or the second law depends on the temperature of transition if the ice temperature at the base of the ice sheet drops below 6 5 C then the second table of parameters loidef_2 are used Ice viscosity for Shallow Shelf Approximation amp diagno_rheol nom du bloc diagno_rheol sf01 0 2 sf03 0 03 pvimin 1 e3 coefficients related to the flow law of grounded ice sf01 viscosity coefficient linear law sf03 viscosity coefficient law n 3 pvimin pvi value for ghost nodes 1 e3 very small relative to standard values 1 e10 This blocks is linked with the previou
50. he following relationship coyrlul Tuz Th m L m er Ts Tr 51 where m corresponds to the melting rate co 3974 Jkg C and cr 2009 Jkg C are the specific heat capacities of water and ice respectively u and Tmz are the speed and GRenoble Ice Shelf and Land Ice model a practical user guide temperature of the ocean mixed layer T is the temperature of the ice ocean interface i e at the base of the ice shelf yr is a coefficient representing the transfer of heat and salt through the boundary layer L 3 35 x 10 Jkg is the latent heat of ice fusion T 25 C is the core temperature of the ice shelf In their ocean model the ice shelf basal temperatures T is calculated using a relationship dependent on salinity not shown In GRISLI we use a vertical ocean temperature profile T from a climate simulation to interpolate the temperature at the base of the ice shelf By doing so the computed vertical temperature profile already account for the variation in temperature and salinity with depth However this does not account for the variation in salinity caused by the melting of the ice shelf itself Our approach here is therefore simplified The input vertical profile file is bi dimensional and should have three columns the depth of each vertical level given in centimetres the index number of each level and the temperature at each level extracted at the chosen location The forcing file is located
51. hy grid and climate forcing basal hydrology and basal melting while the others are located with the other routines in GRISLI SOURCES COMPILATION One the user has defined the modules that will be compiled for the geoplace the routines may be compiled through the Makefile as follow gt cd GRISLI SOURCES gt make Hemin 40 Hemin 40 is the name of the executable that will be produced by the compilation of the various routines contained in the Makefile This executable will be placed in GRISLI bin See the next section for a more detailed explanation of the Makefile NAMELIST SETTINGS The setting of the namelist options can be done after the compilation In the following we provide a detailed description of the various blocks of the namelist Each grid domain in GRISLI has its own namelist file located under GRISLI SOURCES Fichiers parametres In the case of the Northern Hemisphere the namelist file is hemin40_param_list dat Centro Euro Mediterraneo sui Cambiamenti Climatici O7 O Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers Run name and restart procedure amp runpar name of the run parameters block runnam M NHtestl 8 characters at maximum icompteur 0 reprcptr o itracebug 0 num_tracebug 166 comment_run Northern Hemisphere test run runname nom de 1 experienc 8 caracteres icompteur restart in a file 0 gt no l 1 gt yes 2 gt only
52. iations in ocean temperature J Climate 21 11 2558 2572 Hutter K 1983 Theoretical glaciology material science of ice and the mechanics of glaciers and ice sheets Reidel Publishing Company Dordrecht The Netherlands Laske G and Masters G 1997 A Global Digital Map of Sediment Thickness Le Meur E and Huybrechts P 1996 A comparison of differ ent ways of dealing with isostasy examples from modeling the Antarctic ice sheet during the last glacial cycle Annals of Glaciol ogy 23 309 317 MacAyeal D 1989 Large Scale Ice Flow Over a Viscous Basal Sediment Theory and Application to Ice Stream B Antarctica J Geophys Res 94 4071 4087 Peyaud V 2006 Role of the Ice Sheet Dynamics in major cli mate changes PhD thesis Laboratoire de Glaciologie et de G ophysique de l Environnement Universit Grenoble Peyaud V Ritz C and Krinner G 2007 Modelling the Early We ichselian Eurasian Ice Sheets role of ice shelves and influence of ice dammed lakes Clim Past 3 375 386 Pfeffer W T Meier M F and Illangasekare T H 1991 Retention of Greenland runoff by refreezing Implications for projected future sea level change J Geophy Res 96 C12 22117 22124 Reeh N 1991 Parameterization of Melt Rate and Surface Tem perature on the Greenland ice Sheet Polarforschung 59 3 113 128 Ritz C Rommelaere V and Dumas C 2001 Modeling the evolution of Antarctic ic
53. ing_heino use use use use use Basal dragging dragging_vitbil dragging_hwatstream dragging_hwat_cont dragging_hwat_contmaj friction streams Cat Antarctique dragging_hwat_sedim friction with sediments maps use dragging_jorge friction stream Jorge Sediments TAR Elliptique equation eq_elliptique_mod old version remplimat 5 eq_ellip_sgbsv_mod new version July 2008 Calving calving_vincent old version calving_frange new version SSS Mass Conservation equat_adv_diff_2D with advection diffusion use use Basal Melting bmelt_ant_regions For Antarctica with regions bmelt_nor_regions For Northern Hemisphere with regions GRenoble Ice Shelf and Land Ice model a practical user guide use bmelt_nor_depth For Northern Hemisphere with water depth criterion use bmelt_eurasie_regions For Eurasia with regions use bmelt_eurasie_depth For Eurasia with water depth criterion Fake routines for compatibility use fake_heino use fake_eurasi use fake_nor Outputs use out_profile ASCII vertical profiles use out_cptr ASCII restart files use out_hz ASCII variable outputs Time outputs use output_hemin40_mod use output_eurasie20_mod Debugging use printtable end module module_choix Some of those routines are located under GRISLI SOURCES Hemin40 files namely all the routines related to initial topograp
54. is 3000 years Turcotte and Schubert 2002 Centro Euro Mediterraneo sui Cambiamenti Climatici Centro Euro Mediterraneo sui Cambiamenti Climatici O D CMCC Research Papers Ice load pst Lithosphere Asthenosphere Figure 8 Sketch of the ELRA method Figure from Greve and Blatter 2009 GRenoble Ice Shelf and Land Ice model a practical user guide Tutorial setting a Northern Hemisphere simulation In this section the setting of a run over a Northern Hemisphere grid is detailed and explained There are five mains step to run a simulation 1 Selecting some aspects of the physics ND Compilation ao Namelist settings 4 Output variables and frequency 5 Running For each of the grid domain implemented in GRISLI a short name related to the geographical domain is set and reused within the routines The Northern Hemisphere domain implemented in GRISLI has been projected using the Lambert Equal Area projection with a 40 km horizontal resolution This domain covers all the Northern mid latitudes from about 20 N to the North Pole The short name attributed to this domain or geoplace in GRISLI is hemin40 and will be used in the following instruction to refer to the domain PHYSICS OPTIONS In the previous section the main aspects of the physics embedded in GRISLI has been described Some of those aspects present several options that may be set in the namelist or in a module related
55. level value used to perturb the sea level In the example shown above the file is set to perform a steady state simulation for a duration of 20 000 years with no temperature and sea level increments Those two values are nevertheless interpolated during run time at GRISLI time step In this case interpolation result is always 0 rappact corresponds to the precipitation factor y used to correct precipitation for elevation in in the accumulation routine retroac is not used anymore the other options are similar to that set in the clim_forcage namelist block Positive Degree Day parameters ablation amp ablation_ann PDD base Tann et Tjuly Centro Euro Mediterraneo sui Cambiamenti Climatici Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers Cice 0 008 Csnow 0 003 Sigma 5 csi 0 6 Cice and Csnow melting factors for ice and snow sigma variabilite Tday csi proportion of melted water that can refreeze This block sets the value of the parameters for the ablation calculation through the Positive Degree Day semi empirical method See the Numeric and Theory section m Cice is the ice melting coefficient Csnow is the snow melting coefficient Sigma is the amplitude of the daily temperature cycle in C csi is the fraction of melting that can refreeze over the ice sheet Coupling between SIA and SSA mechanics meca_SIA_L1 bloc resol_meca i_re
56. loating ice Figure 6 During the simulation the location of these different regions is determined in every time step The three different regions are characterised by distinct flow behaviour 1 Ice sheets inland ice are areas of grounded ice where ice flow results from deformation In addition basal sliding might occur in areas where the base is at the pressure melting point Here the SIA applies 2 Ice streams are regions of fast flowing ice In reality ice streams are glaciological features with widths of a few kilometres much smaller than our grid size 20 to 40 km In GRISLI the effect of these features is parametrised by treating areas which have the large scale characteristic of ice streams differently i e narrow valleys thick sediment layer saturated by melt water high effective basal pressure In these areas we apply the SSA a different set of momentum equations than for the slow flowing ice sheet following the approach by MacAyeal 1989 3 Ice shelves are large floating ice plates in contrast to ice sheets and streams which are grounded and are fed by inland ice The line where ice begins to float is called the grounding line and once the flow passes this line the ice mass does not contribute to sea level any more Ice shelf nodes must satisfy a flotation criterion and velocities in these grid points are computed from the same set of equations as the velocities in ice stream nodes but assuming a zero basal drag NO
57. lt corresponds to the climate related coefficient to modulate the ice shelves basal melting when meth_Hcoup option is set to 1 in the calving namelist block maxcoefbmelt same as for mincoefbmelt rapbmshelf only use for Antarctica modulate the melting of the ice shelves as well m filforc corresponds to the input climate index file E tempgrad is the mean annual atmospheric lapse rate value C m E tempgrjul is the July atmospheric lapse rate value C m Perturbation of the reference climate forcing GRenoble Ice Shelf and Land Ice model a practical user guide amp clim_pert nom du bloc rappact 0 05 0 05 pour Antartique retroac 0 rapbmshelf 5 mincoefbmelt 0 maxcoefbmelt 2 filforc b140_topol_150 dat tempgrad 0 005 Annual lapse rate C m tempgrjul 0 004 Summer lapse rate C m To perform a steady state simulation this module must be commented out in the module choix In this module GRISLI also offers to perturb the reference climate using some time varying values prescribed in the filforc similarly to the previous namelist block However this forcing file is slightly different and must follow the structure nb_lines 2 0 Ow 20000 0 0 where nb_lines corresponds to the number of lines contained in the file The first column is the time the second corresponds to temperature value used to perturb the reference climate in C and the last column is a sea
58. n the next step will be tbegin dtt until tend is reached testdiag is atime step needed to diagonalise the elliptic matrix and modulated by the mass conservation routine In case of crash the user may change this value to adjust the dynamical time step dt calculated within the mass conservation routine The values is also function of the grid resolution and a typical value for a 40 km grid is about 0 02 Grid topography GRenoble Ice Shelf and Land Ice model a practical user guide amp topo_file filel NH_40km_GRISLI xyz file2 NH_40km_GRISLI xyz filecoord Coordinates_NH_40km_GRISLI xyz filel initial topo file2 reference topo filecoord grid coordinates In this block the user can set the file containing the initial and reference topography on the GRISLI grid fi1le1 contains the initial topography interpolated on the GRISLI grid This can corresponds to a climate model topography horizontally interpolated on the GRISLI grid for example but showing a less resolved topography than the GRISLI topography contained in file2 m file2 contains the reference topography of GRISLI which is a high resolution topography filecoord contains the grid coordinates i j X Y lon lat see Section 4 Those two topographic files are used to correct the climate variables during run time As described in Section 5 because the climate fields might not have been computed at the GRISLI
59. nd visualise GRISLI outputs the user might use the following tools only indicative NCVIEW available at http meteora ucsd edu pierce ncview_home_page html This tool allows for quick visualisation of the Netcdf both during run time and after the simulation This tools requires the installation of Netcaf a NCO NetCDF Climate Operators available at http nco sourceforge net This tool allows for direct operations on the output files such as the extraction of a single variable the merging of several files arithmetic basic operations etc CDO Climate Data Operators available at https code zmaw de projects cdo This tool is similar to NCO but might also offer more options to perform high level arithmetic operations on the files NCL NCAR Command Language available at http www ncl ucar edu This open source software developed and maintained by NCAR is an advanced tool to plot and perform high level mathematical calculations on netcdf files This tools handle the Cartesian native grid which are the grids used by GRISLI which allow the user to draw projected maps GRenoble Ice Shelf and Land Ice model a practical user guide GMT Generic Mapping Tools available at http www soest hawaii edu gmt This software is similar to NCL but does not handle very easily Netcdf files and allows for only limited math ematical operations on the variables However this tools handles perfectly all the Cartesian projectio
60. ng factor of snow Precipitation number of positive degree days Surface mass balance Annual mean surface air temperature Mean July surface air temperature Daily surface air temperature Melting rate at the base of the ice shelf fraction of ice that refreezes Lapse rate for annual mean air temperature Lapse rate for July air temperature Precipitation correction factor Standard deviation of the daily air temperature SS ES C C C m yr m yr m yr m yr m yr mm d C mm d C m yr m yr C C C m yr C km C km Cri C Thermal properties Qi E K Ice deformation Bar T Ey SI Sliding and dragging Er N Tb Basal hydrology K Ve Pw hw Physical constants GRenoble Ice Shelf and Land Ice model Table 2 Notation used in this report deformational heat heat capacity of ice thermal conductivity of ice Coefficient of the Glen flow law Enhancement factor of the flow law Stress tensor Deviatoric stress tensor Strain rate exponent in the Glen flow law Vertically integrated viscosity of the ice Basal dragging coefficient Effective pressure at the base of the ice Basal drag Hydraulic conductivity Sub glacial flow velocity Subglacial water pressure Hydraulic head Density of ice Density of sea water Density of fresh water Gravitational acceleration Ideal gas constant a practical user guide Wm Jmol kg Ww oga1 m1 Pa
61. nge 4 gt option 3 but with test on upstream point ice thickness point meth_Hcoup gt Hcoup depends on surface temperature If 0 then use namelist Hcoup value Centro Euro Mediterraneo sui Cambiamenti Climatici O Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers The calving is parametrised according to if the Lagrangian scheme described in the theory above fulfils some of the tests upon the flux and the mass balance of the ice shelves The calving occurs if the ice shelf front thickness reaches the critical value Hcoup This value can be prescribed directly in the namelist or can be determined based on the climate forcing Hcoup threshold for ice shelf front thickness cut in meters E ifrange type of tests performed on the Lagrangian ice flux Option 1 and option 2 are biased and should not be used meth_Hcoup if 1 Hcoup depends on surface temperature lf 0 then the calving routine uses namelist Hcoup value Deformation law loi de deformation 1 amp loidef_1 exposant_1 3 temp_trans_1 6 5 enhanc_fact_1 3 coef_cold_1 1 660 O cold 1 7 8201 coef_warm_1 2 000 Q_warm_1 9 545 E E E E 16 04 16 04 module deformation_mod_2lois exponent glen transition temperature ttrans enhancement factor sf for temperatures lower than Temp_trans coef_cold Batl for temperatures greater than Temp_trans coef_ warm Bat2
62. ns on sphere and ellipsoid It also handles very easily ASCII files Note that only open source software are listed here and this list is by far not comprehensive other tools such as GNUplot R Ferret or GRADS might also work 2 GRISLI directory structure The GRISLI main directory includes six sub directories Figure 1 GRISLI SOURCES contains all the routines of the model the Makefile and the namelists to be set up by the user before running the model GRISLI INPUT contains all the input files needed by GRISLI to initialise a simulation such as the model grids the climate forcing and other forcing files GRISLI RESULTATS contains all the output files in ASCII and NetCDF format GRISLI Fichier CPTR contains the restart file in ASCII format GRISLI Param contains a copy of the namelist options created during the initialisation phase of the simulation GRISLI bin contains the model executable created during the compilation as well as the log file containing all the run time informations When getting the code for the first time GRISLI Param GRISLI bin GRISLI Fichier CPTR and GRISLI RESULTATS does not exit and have to be created by the user 3 Compiling and running Before compiling the routines the user has to edit the Makefile located in GRISLI SOURCES There are four steps to execute STEP 1 NetCDF libraries setting the path of the NetCDF include files and libraries in the variables NCD
63. ould be located under GRISLI INPUT Forcage and should look like nb_lines 71 21000 0 00000000 0 00000000 0 00000000 20900 3 80087309E 02 3 80087309E 02 0 00000000 20800 5 28201908E 02 5 28201908E 02 0 00000000 GF FI where nb_lines corresponds to the total number of lines contained in the file The first column corresponds to the time in years BP the second column contains the temperature index Tindex the third column contains the precipitation index Pindex which is similar to Tindex in this case and the last column contains the sea level index SLindex here 0 which implies that the sea level value will not be varied during run time The climate index only depends on time and is informally applied to the whole grid domain In this example the values in the file have a time step of 100 years GRISLI time step is 1 year therefore during run time the climate index values are also linearly interpolated on GRISLI time step Note that for transient experiments the time red in the headers of the climate forcing snapshots corresponds to the real time of the model The model is able to account for negative paleotimes therefore the namelist variables has to bet set such that tbegin 21000 and tend 0 m forcage_filel corresponds to the climate snapshot expressed as anomaly relative to the reference climate Climate will be interpolated linearly during run time between this climate and the reference climate E mincoefbme
64. quals surface topog raphy and ice thickness if some ice sheets already exist in the initial topography For example if you consider running GRISLI over present day Northern Hemisphere topography you will have to provide Greenland ice thickness and the bedrock topography below Greenland In Figure 2 we provide the example of an initial topography corresponding to a glaciation period over the Eurasian Centro Euro Mediterraneo sui Cambiamenti Climatici CO Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers domain described in the previous section For each of the grid point the three topography variables have to be provided in meters and written in the topography file located under GRISLI INPUT EURASIE20 as follows Eurasia present day topo interpolated by grid maker 56700 3 210 270 20 0000 S H BO 240 116241 0 00000000 240 116241 146 568726 0 00000000 146 568726 126 997734 0 00000000 126 997734 113 883545 0 00000000 113 883545 S corresponds to the surface topography H is the ice thickness and BO is the bedrock topography The values are ordered exactly as in the coordinate file described in the previous section The three first lines of the file are mandatory and are red by GRISLI during the initialisation of the grid In the second line some grid informations are indicated total number of points X Y the number of columns contained in the file number of X points number of Y points horizont
65. s less dense than the sea water and flows toward to surface following the base of the ice shelves This fresh water mass is deviated by the Coriolis strength and exchange of mass temperature and salinity occurs with the surrounding sea water This fresh water mass mixes up completely when its density reaches that of the sea water Basal melting decreases gradually from the grounded line toward the front of the ice shelves and some refreezing can occur in the central part At the front the basal melting is generally high because this part of the ice shelf is located far away from the ice sheets and is therefore subject to warmer oceanic currents While several parametrisations of basal melting under the ice shelves have been developed in the ocean general circulation models most of the ice sheet models are not coupled with ocean therefore it is difficult to properly calculate the basal melting accounting for ocean temperature changes In GRISLI basal melting follows the work of Dumas 2002 who calculates the basal melting in function of surface temperature and an ocean depth criterion A distinction is made between the grounded line bmgrz and the rest of the ice shelf bmshe11 1 according to the depth criterion set in the namelist two different values of basal melting can be attributed one above this limit bme1t1 and one below bme1t2 During the initialisation phase bmshel f bmeltl 39 bmgrz bmelt1 40 or bmelt2 if appropriate
66. s one and is part of the Glen flow law calculation see Numeric and Theory It calculates the thermo mechanically coupled integrated ice viscosities between the grounded ice SIA and the floating ice SSA sf01 and sf03 correspond to the enhancement factors of the floating part Basal Melting under ice shelves amp bmelt_nor_depth module bmelt_nor_depth m bmelt_type bmelt1 0 bmelt2 0 depthlimit 2000 ocn_profile TEMP_vert_Chukchi dat ocn_map Tocn_B140_topol_artic_20km_noice_GRISLI xyz bmelt_type 1 gt prescribed basal melting depending on depthlimit Centro Euro Mediterraneo sui Cambiamenti Climatici O1 Centro Euro Mediterraneo sui Cambiamenti Climatici CMCC Research Papers bmelt_type 2 gt basal melting calculated using an ocean temperature vertical profile bmelt_type 3 gt as for 2 but using bidimensional ocean vertical temperature bmelt1 basal melting above depthlimit bmelt2 basal melting below depthlimit depthlimit depth threshold to precribe basal melting ocn_profile ocean temperature vertical profile ocn_map vertical ocean temperature map sbmelt_nor_regl nom du bloc lbmelt_Ross 0 5 bmgrz_Ross 0 2 lbmelt_FRis 2 4 bmgrz_FRis 5 6 In the namelist two blocks are associated with ice shelves basal melting namely bmelt_nor_depth and bmelt_nor_reg1 The choice of the basal melting module is operated in the module_choi
67. search Papers false true annual type_precip filel reference climate annual true if climate file contains Tann Tjja type_precip true if liquid precip false snow ij amp clim_pert nom du bloc rappact retroac rapbmshelf mincoefbmelt maxcoefbmelt filforc tempgrad tempgrjul 05 0 05 pour Antartique NOudoo b140_topol_150 dat 0 005 Annual lapse rate C m 0 004 Summer lapse rate C m sablation_ann PDD base Tann et Tjuly Cice Csnow Sigma csi ouoo Cice and Csnow melting factors for ice and snow sigma variabilite Tday csi proportion of melted water that can refreeze meca_SIA_L1 bloc resol_meca i_resolmeca 2 i_resolmeca type d association entre SIA et L1 i_resolmeca 0 chacun dans sa zone i_resolmeca 1 dans les zones stream addition si uxdef gt uxLl MIS11 Cairns i_resolmeca 2 addition systematique dans les zones stream 8 Output variables and output frequency The output variables that have to be written during run time in the output NetCDF and ASCII files are listed in the file LISTE VAR NETCDF dat and LISTE VAR HZ dat respectively Those files are located under GRISLI SOURCES Fichiers parametres The variables are listed according to their dimension and the user may decide the frequency of output of those variables by switching on some flags For example S name of the variable in GRISLI index number of the v
68. solmeca 2 i_resolmeca type d association entre SIA et Ll i_resolmeca 0 no overlapping between the two mechanics ji_resolmeca 1 in the stream areas addition if uxdef gt uxLl i_resolmeca 2 systematic addition in stream areas This option is usually set on 2 GRenoble Ice Shelf and Land Ice model a practical user guide Tutorial how to implement a new grid domain To implement a new grid domain in GRISLI there are several steps to follow 1 create a new input repository in GRISLI INPUT to put the grid the geothermal heat flux the sediment map and the climate forcing interpolated on the new grid For example GRISLI INPUT NEWGRID 2 create a new repository in GRISLI SOURCES copying the routines located under one of the pre existing implemented domain into the new grid repository For example GRISLI SOURCES NEWGRID 3 editthe various routines to give the dimension of the grid and point to the right input boundary and initial conditions files 4 edit the Makefile located under GRISLI SOURCES to add the new domain compile the corresponding routines and create the executable In this guide the GRISLI repository contains two different grid domains one covering the entire Northern Hemisphere from about 16 N and one covering Northern Eurasia The grid and climate forcing related to those grids are located under GRISLI INPUT HEMIN40 and GRISLI INPUT EURASIE20 The routines specific to those two domains ar
69. surface elevation at which the climate forcing was computed or interpolated to in case of reanalysis presents some discrepancies with the finer surface topography of the ice sheet model Therefore a correction of the climate variables is necessary before giving it to the ice sheet model for simulation Tocem EMB ES GCM ism eW TGCM TISM GCM Figure 3 Example of downscaling from a coarse climate model GCM onto a finer resolution ice sheet model ISM Initial GCM topography Sac m is less resolved than initial ice sheet topography Sz sm and in this example Sac y is lower than S7 51 To correct the simulated GCM temperatures Tac m the downscaling is performed by means of an atmospheric lapse rate A Simulated GCM precipitations Paga m can be corrected following an exponential relationship accounting for the temperature correction due to elevation changes This procedure is the standard procedure applied to run an ISM However more sophisticated methods could be use for example to conserve the mass during the interpolation Downscaling procedure 1 The temperature field from the CESM simulation is horizontally interpolated on the higher resolution Eurasian ice sheet model ISM grid Tecm 1 x 1 gt Tocm on 1sm 20km x 20 km 1 2 Interpolation of the surface elevation field from the coarse resolution CESM grid onto the ice sheet model grid introduces small differences in elevation Temper
70. surface temperature If 0 then use namelist Hcoup value 1 1 loi de deformation 1 module deformation_mod_2lois amp loidef_1 w exposant_1 temp_trans_1 6 5 enhanc_fact_1 3 coef_cold_1 1 660E 16 Q_cold_1 7 820E 04 coef_warm_1 2 000E 16 Q_warm_1 9 545E 04 exposant glen temperature de transition ttrans enhancement factor sf pour les temperatures inf a Temp_trans coef_cold Bat1 et Q cold Q1 pour les temperatures sup a Temp_trans f coef_warm Bat2 et Q warm 02 1 1 loi de deformation 2 module deformation_mod_2lois amp loidef_2 exposant_2 temp_trans_2 10 I Ha GRenoble Ice Shelf and Land Ice model a practical user guide enhanc_fact_2 10 coef_cold_2 8 313E 08 Q_cold_2 4 000E 04 coef_warm_2 8 313E 08 Q_warm_2 6 000E 04 enhancement factor sf pour les temperatures inf a Temp_trans coef_cold Batl et Q_cold 01 pour les temperatures sup a Temp_trans exposant glen temperature de transition ttrans 1 1 y coef_warm Bat2 et Q warm 02 1 amp diagno_rheol nom du bloc diagno_rheol sf01 0 2 sf03 0 03 pvimin 1 e3 coefficients par rapport a la loi glace posee sf01 coefficient viscosite loi lineaire sf03 coefficient viscosite loi n 3 pvimin valeur de pvi pour les noeuds fictifs 1 e3 tres petit par rapport aux valeurs standards 1 e10 amp bmelt_nor_depth module bmelt_nor_
71. t This process is computed at each time step and when the thickness of the ice shelf front is lower than the calving criterion GRISLI tests the status of the front The ice shelf thickness variations are computed by means of a Lagrangian scheme Ou Ov bm bmelt mo dH dt 34 where H is the ice thickness u and v are horizontal velocities bm corresponds to the surface mass balance and bmelt corresponds to the fusion at the base of the ice shelf prescribed in the namelist This equation is applied to the upstream node of the pixel tested for calving The upstream node has an ice thickness Hupstream greater than the thickness criterion Hcoup The time used by the ice in the upstream node to flow to the neighbour ice shelf front point is given by Ox gt 35 gt 35 Therefore the ice thickness variation during this time 7 is dh Tp T 36 This quantity is generally negative because the upstream point lose some mass in favour of the ice shelf front node Using this quantity the thickness at the front H front can be deduced dh H front Hupstream Tf dt 37 If H front is lower than H oup calving occurs in the node In the case that no calving occurs the condition to maintain and develop the ice shelf front is given by dH dt lupstream gt Hcoup Hupstream 38 de As reported by Peyaud 2006 a physical problem arises when ice shelves develop from coastal grounded ice In the coastal
72. t the sediment layer has to be thick ha gt 150 m However this threshold is ad hoc and no direct evidence exists to constrain it better see Discussion Second it has to be saturated by meltwater This saturation is determined by the hydraulic head hw which has to exceed 250 m the hydraulic head is calculated with a Darcy flow law see section Third the effective pressure N has to be lower than 3 5 107 Pa The grid point is located on an ice rise with low ice thickness and low surface slope lt has to hold that pg HVS lt 7 10 Pa BASAL HYDROLOGY One criterion for determining the position of ice streams is the basal water head h In the following we describe how basal hydrology is implemented in GRISLI The basal hydrology is characterised by different hydrological features ranging from drainage through a thin film of water lt 1 mm to flow through a network of cavities and rivers Also percolation of basal water into the porous sub glacial sediment plays an important role Nonetheless due to a lack of observations we still do not fully understand the hydrological processes occurring below ice sheets Thus the representation of basal hydrology in an ice sheet model is one of the most arguable components but indeed one of the most important processes for the dynamics of the whole ice sheet The basal hydrology in GRISLI was implemented by Peyaud 2006 and is based on a Darcy type flow law The drainage of sub glacial wa
73. ter depends on the pressure imposed onto the base This pressure can be expressed as a potential and the flow of melt water follows its gradient The potential can be written as py Pfw9B pgH Bo 28 Pw Stands for the sub glacial water pressure The second term represents the influence of the altitude of the base while the third term represents the pressure resulting from the weight of the ice Po is an arbitrary constant Then the gradient of can be expressed as VO Vow V PfugB pgH Vow pg ors 1 VB VS 29 Note that the magnitude of the surface gradient VS is around 10 times larger than that of the bedrock Darcy s law states that flow through a porous medium is proportional to a hydraulic gradient which corresponds to the variation of the hydraulic head hw This leads to an equation for the flow velocity and with expressing the hydraulic head as h pywg we obtain K Pwg Ve KVhy VO 30 with the proportionality constant K the hydraulic conductivity This parameter depends on the property of the soil e g clay lt 107 m s till 1071 to 1075 m s fine sand 1077 to 1075 m s In GRISLI K is set in the following way If the base is cold with the ice temperature below the GRenoble Ice Shelf and Land Ice model a practical user guide pressure melting point we make the assumption that the ground is impermeable and thus K is set to zero If the base is warm with the ice temperature a
74. to simulate both grounded and floating ice The grounded part uses the Shallow Ice Approximation SIA Hutter 1983 whereas ice shelves and ice streams are simulated following the Shallow Shelf Approximation SSA MacAyeal 1989 The ice shelf formulation in GRISLI allows for a more realistic calculation of the ice sheet growth and particularly of the advance of ice onto the shallow continental shelves in both hemispheres high latitudes GRISLI has been validated over Antarctica Ritz et al 2001 and successfully applied to study the inception of ice sheets during the Early Weichselian period Peyaud et al 2007 Alvarez Solas et al 2011 In the following document we give a detailed description of the model and some indication to set up and run simulations This guide is by far non exhaustive and is related to a version of GRISLI that might not be up to date The purpose of this guide is purely descriptive and the implementation performed at CMCC are indicated in the text Finally GRISLI is not open source therefore CMCC is not authorised to diffuse the code The code is available on request to Catherine Ritz LGGE Grenoble 0 Centro Euro Mediterraneo sui Cambiamenti Climatici CD Centro Euro Mediterraneo sui Cambiamenti Climatici D CMCC Research Papers User Guide 1 Before installing GRISLI is a suite of Fortran 90 routines that need three distinct numerical tools to compile and run properly Before running the model th
75. to the grid domain and in which some of the parametrisation can be selected In the case of the Northern Hemisphere domain the file is called module_choix hemin40 0 4 f90 and is located under GRISLI SOURCES Hemin40_files This module allows to choose some physical and practical aspect of the model MODULE MODULE_CHOIX Reading initial topography use lect_topo_anteis For Antarctica 40 km use lect_topo_hemin40 For Northern Hemisphere 40 km use lect_topo_eurasie For Eurasia toute resolutions Reference Climate typ use lect_clim_act_anteis use lect_clim_act_hemin40 For Northern Hemispher use lect_clim_act_eura20 For Eurasia Climate forcing typ Centro Euro Mediterraneo sui Cambiamenti Climatici GO Centro Euro Mediterraneo sui Cambiamenti Climatici O D CMCC Research Papers use use use climat_perturb_mod climat_forcage_mod climat_profil_mod use climat_regions_delta See Proglacial lakes use lakes_prescribed_ mod no_lakes use Isostasy isostasie mod module accounting for isostaty NOISOSTASIE MOD no isostasy use Ice deformation law deformation_mod deformation law use use Ice thermal properties prop_thermiques_ice prop_therm_ice_heino use Basal water hydrology eau_basale use use use Basal sliding module sliding_vitbal sliding_Bindschadler sliding_dragg
76. x The first block can be use for any grid domain except Antarctica while the second block is related to Antarctica grid domain exclusively In the case of Antarctica bmelt_nor_reg1 attributes a melting value to some specific ice shelves grounded line and central area Dumas 2002 In the block bmelt_nor_depth on the contrary three methods two out of three implemented at CMCC to prescribe basal melting can be chosen E bmelt_type define the way of accounting for ocean heat fluxes bmelt_type 1 corre sponds to the standard method to prescribe basal melting depending on a depth limit The depthlimit option has to be filled with an ocean depth in meters and two basal melting values bmelt1 applied for depths above the prescribed limit and bme1t2 applied for depth below the limit can be set in m yr M bmelt_type 2 is the basal melting calculated using the quadratic parametrisation of Holland et al 2008 in which a vertical ocean temperature profile is prescribed In their paper Holland et al 2008 review the different parametrisations used in the Ocean General Circulation Model to calculate the ice melting under the ice shelves in Antarctica and its interaction with the surrounding ocean waters in terms of temperature and salinity Comparing observations with the various peer reviewed parametrisations they show that the relationship between basal melt rates and vertical ocean temperature below the shelves is quadratic They propose t

GRenoble Ice-Shelf and Land-Ice model: a practical user guide

Contents

Download Pdf Manuals

Related Search

Related Contents