Glacier table

The glacier table allows for a sub-grid representation of the glacier processes, by creating a grid at a finer resolution than the model grid. The glacier table is only required for the SPHY version higher 3.0 and higher. The glacier table contains eight elements (Table 10, Table 11).

Table 10: Elements of the glacier table

Name	Description
U_ID	Unique ID that is set at smaller spatial scale to account for sub grid heterogeneity in the glacier cover
MOD_ID	Individual ID of each cell of the SPHY model
GLAC_ID	ID of the individual glaciers in the model area, from the RGI inventory, from glacier ID 1 to X (with X being the maximum number of glacier)
MOD_H	Elevation of model ID cell
MOD_GLAC	Elevation of subgrid cell
DEBRIS	A flag (0 or 1) to identify if the cell is debris-covered
FRAC_GLAC	The fraction of the cell that is covered by glacier
ICE_DEPTH	Ice thickness of the cell

Table 11: Example of Glacier table

U_ID	MOD_ID	GLAC_ID	MOD_H	GLAC_H	DEBRIS	FRAC_GLAC	ICE_DEPTH
1	61122	89	3946	3946	1	0.002	0
2	61123	89	3928	3928	1	0.368	5.28
3	61124	89	3897	3897	1	0.189	7.04
4	61937	89	3950	3950	1	0.446	17.31
...	...	...	...	...	...	...	...

The suggested method to process these input can be summarized in 4 steps:

1. Create model vector grid (fishnet) (MOD_ID)

2. Intersect with the glacier outlines (GLAC_ID)

3. Create subgrid variability grid (U_ID)

4. Calculate elevation, glacier thickness, debris cover, and glacier fraction (MOD_H, GLAC_H, DEBRIS, FRAC_GLAC, ICE_DEPTH)

Step 1: Create model vector grid (MOD_ID)

Open a model input raster file with the spatial extent and grid size in QGIS. This could be the clone.map, dem.map, latitude.map, or any raster model grid input. In the Vector Creation toolbox, select the “Create grid” tool. For the “Grid extent”, select the option “Calculate from layer”, and select the model raster. For “Horizontal spacing” and “Vertical spacing”, set the grid size of your model. Set name to “MOD_ID”.

Figure 42: Creating Model_ID grid. Set the Horizontal and Vertical spacing to the model grid resolution.

This new vector layer contains polygons for each grid of your model domain, with value starting at 1 in the upper left corner. In the attribute table of your “MOD_ID” layer (right click on layer, select “Open attribute table” and deleted the attributes “left”, ‘”top”, ‘”right”, “bottom” to keep only the column “id” using the “Delete field” in the attribute table. Rename this attribute to “MOD_ID” using the tool “Rename field” in the Vector table toolbox. This is now your MOD_ID vector layer.

Step 2: Create GLAC_ID, MOD_H and glac_id.map

Import the vector layer of glacier outline in the model area. This could come from manually digitized glacier areas or from an inventory like RGI7 (RGI 7.0 Consortium, 2023). Clip the glacier outline vector to the model area using the Vector Overlay toolbox “Clip” tool and reproject to the desired projection to match your model projection (toolbox Vector general, “reproject layer” tool). Open the attribute table of the glacier outline layer and create a new attribute using the field calculator (Figure 43).

Using the GDAL toolbox, use the “Rasterize (vector to raster) to transform the GLAC_ID vector layer to a raster layer GLAC_ID” (Figure 44), followed by transforming to a .map file by using the PCRaster toolbox with ther tool “Convert to PCRaster format”, with format " NOMINAL” and named glac_id.map. This is an input for the glacier module in SPHY.

Figure 43: Create GLAC_ID attribute

Figure 44: Rasterize glacier outline to obtain GLAC_ID map layer.

You now need to intersect the MOD_ID grid with the glacier outline to obtain a grid that covers only the glacier area. To do so, intersect the glacier outline vector layer with the MODEL_ID grid (Vector overlay toolbox, “intersection” tool) with input layer “MOD_ID” and overlay layer “Glacier Outline”. The resulting layer (GLAC_ID) has two attributes: MOD_ID, from the original grid of the model domain, and “GLAC_ID”, the value of the individual glaciers (Figure 45).

Figure 45: The GLAC_ID layer (in pink) with two attributes, GLAC_ID and MOD_ID as shown in the attribute table for a cell on the glacier table for a cell on the glacier toe highlighted in yellow, with MOD ID 839 and GLAC_ID 1. The orange grid in the background correspond to the model ID grid. At this stage, the cell size of the GLAC_ID and MOD_ID are the same size.

The elevation of the MOD_ID cells over the glacier area, which corresponds to the column “MOD_H” in the glacier table, can now be calculated. Import the DEM for your model area. Then, open the tool “Zonal Statistics” in the toolbox Raster analysis (Figure 46). For “Statistics to calculate”, select “mean”. Using the tool “rename field”, rename the resulting field to “MOD_H”. The resulting raster has three attributes: MOD_ID, GLAC_ID, MOD_H.

Figure 46: Using the "Zonal statistics" tool to calculate the elevation of the glacier grid.

Step 3: Create sub-grid variability (U_ID)

In this step, the smaller grid to account to calculate glacier processes at a higher resolution is created. Also, in some cases, some grid cells have 2 different glaciers within them, so we have to create a smaller gid system to account for this sub-grid variability (Figure 47). To create the U_ID (the smaller scale grid), follow the same steps as to create the MOD_ID (step 1), but select a smaller spatial resolution as needed. Create a new attribute, named U_ID, to start the values at 1 (same as Figure 43 to obtain GLAC_ID, using the rows as index). Keep only the attribute U_ID in the vector layer. Intersect the U_ID with the MOD_ID vector layer to keep only the smaller grid in the glacier area (using the “Intersect” tool). The resulting vector layer U_ID has four attributes: MOD_ID, GLAC_ID, MOD_H, and U_ID (Figure 48).

Figure 47: The need for a smaller grid to account for the multiple glacier in one cell model

Figure 48: The U_ID layer showing the smaller scale grid cell, with the four attributes (MOD_ID, GLAC_ID, MOD_H, and U_ID). The orange background correspond to the larger MOD_ID grid.

Step 4: Calculate the elevation (U_MOD), glacier thickness (ICE_DEPTH), debris cover (DEBRIS) and glacier fraction (GLAC_FRAC):

In this step, we use the finer grid U_ID to calculate the values for each grid cell for the various dataset: elevation, the presence of debris, the ice thickness and the fraction glacier cover.

For the H_GLAC (elevation of the glacier), follow the same step as for the MOD_H calculation with the Zonal Statistics tool (Figure 46) using the original DEM, or an updated, finer resolution DEM is needed. Once again, you can rename the field to ICE_DEPTH and remove any attribute that are superfluous.

For glacier thickness (ICE_DEPTH), first load your ice thickness data into QGIS, then use the Zonal Statistic tool. If the ice thickness dataset is not continuous, or comes with multiple files, first merge into a single raster using toolbox SAGA Next Gen (needs to be downloaded as a plugin) raster tools with tool Mosaicking. Set “Overlapping areas” to “maximum”, and create a text file list with all the files to be used in the mosaic. Then using the “Zonal Statistics” tools to calculate the ice thickness for each U_ID cell. Once again, you can rename the field to ICE_DEPTH and remove any attribute that are superfluous.

Some ice thickness might show as zero due to some misalignment of the product of a missing glacier in the ice thickness database. For these cases, you can set a fixed value as the glacier thickness. To do so, in the attribute table, open the field calculator, and for the ICE_DEPTH field, set the expression: if("ICE_DEPTH" is NULL, 25, " ICE_DEPTH") (Figure 49). Repeat the procedure for ice_depth values that are 0.

Figure 49: Replacing "NULL" values with ice thickness of 25m, using the layer "_MEAN", THE OUTPUT OF THE Zonal Statistics tool.

For the debris-cover, import the debris-cover file, and if in vector format, rasterized using the GDAL toolbox, with tool “Rasterize (vector to raster), with the same resolution as the U_ID. Set the burn in value to 1. Using the “Zonal statistics” tool, calculate the presence of debris for each pixel of the U_ID vector. Using the field calculator, change the “NULL” values to “0” (if("_mean" is NULL, 0, "_mean"), as in Figure 49 DEBRIS should be either 0 or 1.

For the fraction of glacier cover in a given cell (GLAC_FRAC), first calculate the area of each polygon with the Vector Geometry toolbox, with tool “add geometry attribute”. This adds area and perimeter to the attribute table. Based on the cell size, you can now calculate the ratio of the cell covered by the glacier. For example, in this case, if the cell is 25 x 25m, the maximum area is 625. Using the field calculator, you can divide the “area” attribute by the maximum area of each cell, giving you the glacier fraction. GLAC_FRAC should be higher than 0 and maximum of 1.

You should now have a vector layer with eight attributes corresponding to the the 8 columns needed in the glacier table, with the spatial resolution of the smaller grid capturing the glacier processes at higher resolution. Once all the calculation are done and checked, export file as .csv, organize column titles and order, and this gives you the glacier table needed as a SPHY input (Figure 50).

Figure 50: Exporting the vector layer as glac_table.csv