From ce009cc68f0b9a6bfe4cf245e58a508e3b292780 Mon Sep 17 00:00:00 2001 From: Chloe Hancock Date: Tue, 31 Mar 2026 16:47:55 +0100 Subject: [PATCH 1/5] Add runoff sensitivty notebook --- _toc.yml | 1 + notebooks/tutorials/runoff_sensitivity.ipynb | 1042 ++++++++++++++++++ notebooks/welcome.ipynb | 7 +- 3 files changed, 1049 insertions(+), 1 deletion(-) create mode 100644 notebooks/tutorials/runoff_sensitivity.ipynb diff --git a/_toc.yml b/_toc.yml index 174698a8..712a20ee 100644 --- a/_toc.yml +++ b/_toc.yml @@ -30,6 +30,7 @@ chapters: - file: book/hydro sections: - file: notebooks/tutorials/hydrological_output + - file: notebooks/tutorials/runoff_sensitivity - url: https://oggm.org/oggm-edu-notebooks/oggm-edu/glacier_water_resources.html title: Glaciers as water resources (OGGM-Edu part 1 - idealized) - url: https://oggm.org/oggm-edu-notebooks/oggm-edu/glacier_water_resources_projections.html diff --git a/notebooks/tutorials/runoff_sensitivity.ipynb b/notebooks/tutorials/runoff_sensitivity.ipynb new file mode 100644 index 00000000..0184422f --- /dev/null +++ b/notebooks/tutorials/runoff_sensitivity.ipynb @@ -0,0 +1,1042 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d1abcbfd", + "metadata": {}, + "source": [ + "## Set Up\n", + "\n", + "First install required packages to run this tutorial. \n", + "\n", + "We will be using level 4 data to allow for the dynamical spinup and the calibration capabilities." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b6e8193", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "from oggm import cfg, utils, workflow, tasks\n", + "from oggm.core import massbalance\n", + "from oggm.core.massbalance import mb_calibration_from_scalar_mb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85804d49", + "metadata": {}, + "outputs": [], + "source": [ + "cfg.initialize(logging_level='WARNING')\n", + "cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-runoff', reset=True)\n", + "cfg.PARAMS['store_model_geometry'] = True" + ] + }, + { + "cell_type": "markdown", + "id": "9d2efb09", + "metadata": {}, + "source": [ + "We start from a well known glacier in the Austrian Alps, Hintereisferner. But you can choose any other glacier, e.g. from [this list](https://github.com/OGGM/oggm-sample-data/blob/master/wgms/rgi_wgms_links_20220112.csv)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "402e44d9", + "metadata": {}, + "outputs": [], + "source": [ + "# Hintereisferner\n", + "rgi_id = 'RGI60-11.00897'\n", + "\n", + "# We pick the elevation-bands glaciers because they run a bit faster - but they create more step changes in the area outputs\n", + "from oggm import DEFAULT_BASE_URL\n", + "gdir_hef = workflow.init_glacier_directories([rgi_id], from_prepro_level=4, prepro_border=80, prepro_base_url=DEFAULT_BASE_URL)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "78de2a38", + "metadata": {}, + "source": [ + "## Generating the Hydrological and Glaciological outputs" + ] + }, + { + "cell_type": "markdown", + "id": "041261fe", + "metadata": {}, + "source": [ + "The pre-processed directories don't automatically have hydrological outputs, so lets run the `run_with_hydro` task below to calculate these! This requires a dynamical spinup, let's choose the time period 2000-2020." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36a47e5d", + "metadata": {}, + "outputs": [], + "source": [ + "# file identifier where the model output is saved\n", + "file_id = '_default'\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # which climate scenario? See following notebook for other examples\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix='_default')) as ds:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds = ds.isel(time=slice(0, -1)).load()" + ] + }, + { + "cell_type": "markdown", + "id": "ac645b31", + "metadata": {}, + "source": [ + "Now that we have the hydrological and glaciological outputs for our desired time period, lets calculate the total runoff for the glacier. There are four key variables that form the total runoff when summed together, and the process to calculate the runoff can be seen below.\n", + "\n", + "The total runoff from the glacier is:\n", + "`runoff` = `melt_off_glacier` + `melt_on_glacier` + `liq_prcp_off_glacier` + `liq_prcp_on_glacier`\n" + ] + }, + { + "cell_type": "markdown", + "id": "1d20120d", + "metadata": {}, + "source": [ + "The figure below shows the runoff for the default calibrated parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d27a8a7d", + "metadata": {}, + "outputs": [], + "source": [ + "sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]\n", + "df_annual = ds[sel_vars].to_dataframe()\n", + "\n", + "# These summed variabels give the total runoff from the glacier\n", + "runoff_vars = ['melt_off_glacier', 'melt_on_glacier','liq_prcp_off_glacier', 'liq_prcp_on_glacier']\n", + "\n", + "# Convert them to megatonnes (instead of kg)\n", + "df_runoff = df_annual[runoff_vars] * 1e-9\n", + "fig, ax = plt.subplots(figsize=(10, 3.5), sharex=True)\n", + "df_runoff.sum(axis=1).plot(ax=ax);\n", + "plt.ylabel('Runoff'); plt.xlabel('Years'); plt.title(f'Total annual runoff for {rgi_id}');\n" + ] + }, + { + "cell_type": "markdown", + "id": "92501cfa", + "metadata": {}, + "source": [ + "This was done with the following parameter set:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a217ac7", + "metadata": {}, + "outputs": [], + "source": [ + "gdir_hef.read_json('mb_calib')" + ] + }, + { + "cell_type": "markdown", + "id": "33a812b6", + "metadata": {}, + "source": [ + "## Parameter Calibration" + ] + }, + { + "cell_type": "markdown", + "id": "0d96ff15", + "metadata": {}, + "source": [ + "The mass-balance parameters at Level 4 are pre-calibrated automatically by the OGGM and stored in the dataframe. However there are flexible mass-balance calibration schemes in the OGGM.\n", + "\n", + "In this tutorial, our aim is to calibrate the mass-balance parameters and investigate the impact that this calibration has on the run-off. In the OGGM we can currently (at v1.6.3) calibrate the `melt_f`, `prcp_fac` and `temp_bias` parameters. We will calibrate these parameters using the **scalar mass balance** method, and investigate the relationship with the runoff components and the total runoff output." + ] + }, + { + "cell_type": "markdown", + "id": "ac6c6eaa", + "metadata": {}, + "source": [ + "### **Calibration:** Scalar Mass Balance Calibration" + ] + }, + { + "cell_type": "markdown", + "id": "d5f9d1ac", + "metadata": {}, + "source": [ + "Firstly comparing the mass balances below for in-situ observations and the default calibrated parameters. To understand the current behaviour of the mass balance output.\n" + ] + }, + { + "cell_type": "markdown", + "id": "fe1bd5be", + "metadata": {}, + "source": [ + "Now let's experiment by calibrating the parameters: these parameters are the melt factor (`melt_f`), the precipitation factor (`prcp_fac`) and the temperature bias (`temp_bias`). We will calibrate each parameter and permutations of the calibrations:\n", + "1. Just calibrating `melt_f`\n", + "2. Just calibrating `prcp_fac`\n", + "3. Just calibrating `temp_bias`\n", + "4. Calibrating `melt_f` and `prcp_fac`\n", + "5. Calibrating `prcp_fac` and `temp_bias`\n", + "6. Calibrating `melt_f` and `temp_bias`\n", + "7. Calibrating all parameters; `melt_f`, `prcp_fac` and `temp_bias`\n", + "\n", + "The next calibration section will be relatively repetative, but take note of the parameters being calibrated and the inputs being used in the calibration process.\n", + "\n", + "We will start the calibration by just calibrating the model the `melt_f` paramater." + ] + }, + { + "cell_type": "markdown", + "id": "33360649", + "metadata": {}, + "source": [ + "
\n", + " \n", + " Note: In the OGGM Scalar Calibration method, even when all parameters are nominated for calibration, they may not reach the calibration stage. The OGGM works by calibrating one parameter at a time, only moving onto the next when the previous parameter calibration cannot reach the desired target. In the tutorial below, we will see the result of this in action with usually only the first parameter being calibrated, even if all parameters are nominated for calibration.\n", + " \n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "f47b6995", + "metadata": {}, + "source": [ + "### Single Parameter Calibration" + ] + }, + { + "cell_type": "markdown", + "id": "d28b4228", + "metadata": {}, + "source": [ + "Fetch the reference mass balance we want to calibrate to from the geodetic mass-balance data, this is the data from Hugonnet et al 2021:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e675125", + "metadata": {}, + "outputs": [], + "source": [ + "ref_mb_df = utils.get_geodetic_mb_dataframe().loc[gdir_hef.rgi_id]\n", + "ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == cfg.PARAMS['geodetic_mb_period']].iloc[0]\n", + "# dmdtda: in meters water-equivalent per year -> we convert to kg m-2 yr-1\n", + "ref_mb = ref_mb_df['dmdtda'] * 1000\n", + "ref_mb" + ] + }, + { + "cell_type": "markdown", + "id": "3ec4b550", + "metadata": {}, + "source": [ + "Now we will calibrate with the scalar mass balance and use its default calibration setting. This is just calibrating the `melt_f` parameter:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b5dab26", + "metadata": {}, + "outputs": [], + "source": [ + "# Just calibrate the melt_f parameter\n", + "calib_param_melt_f = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " overwrite_gdir=True)\n", + "\n", + "# Creating a mass balance data frame to store the results\n", + "mbdf= pd.DataFrame(index = np.arange(2000,2020,1))\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['melt_f_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "calib_param_melt_f" + ] + }, + { + "cell_type": "markdown", + "id": "b2f0fcc5", + "metadata": {}, + "source": [ + "We can see the parameters above, we will use the run_with_hydro task to gather the hydrological outputs. We will now calculate the runoff from the calabrated parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97847d11", + "metadata": {}, + "outputs": [], + "source": [ + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "file_id = '_melt_f'\n", + "\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_melt_f:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_melt_f = ds_melt_f.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated melt_f parameter\n", + "sel_vars = [v for v in ds_melt_f.variables if 'month_2d' not in ds_melt_f[v].dims]\n", + "df_annual_melt_f = ds_melt_f[sel_vars].to_dataframe()\n" + ] + }, + { + "cell_type": "markdown", + "id": "2f0173ff", + "metadata": {}, + "source": [ + "Now calibrating the precipitation factor variable, `prcp_fac`. Here we fix the `melt_f` and the `temp_bias` and then calibrate the precipitation factor accordingly, this will ensure that the reference average mass balance will be matched." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4182717", + "metadata": {}, + "outputs": [], + "source": [ + "# Just calibrate the prcp_fac parameter\n", + "calib_param_prcp_fac = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " calibrate_param1='prcp_fac',\n", + " overwrite_gdir=True)\n", + "\n", + "file_id = '_prcp_fac'\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['prcp_fac_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_prcp_fac:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_prcp_fac = ds_prcp_fac.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated melt_f parameter\n", + "sel_vars = [v for v in ds_prcp_fac.variables if 'month_2d' not in ds_prcp_fac[v].dims]\n", + "df_annual_prcp_fac = ds_prcp_fac[sel_vars].to_dataframe()\n", + "\n", + "calib_param_prcp_fac" + ] + }, + { + "cell_type": "markdown", + "id": "86551587", + "metadata": {}, + "source": [ + "Note here that now the `melt_f` is lower than before, and the `prcp_fac` is now lowered in calibration so the average reference mass-balance can be matched." + ] + }, + { + "cell_type": "markdown", + "id": "30085bd7", + "metadata": {}, + "source": [ + "Now calibrating the temperature bias parameter, `temp_bias`. Here we now fix the `melt_f` and `prcp_fac` to the default values, and calibrate the `temp_bias` to match the average reference mass-balance:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0279da19", + "metadata": {}, + "outputs": [], + "source": [ + "# Just calibrate the temp_bias parameters\n", + "temp_bias_calib = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " calibrate_param1='temp_bias',\n", + " overwrite_gdir=True)\n", + "\n", + "file_id = '_temp_bias'\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['temp_bias_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_temp_bias:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_temp_bias = ds_temp_bias.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated melt_f parameter\n", + "sel_vars = [v for v in ds_temp_bias.variables if 'month_2d' not in ds_temp_bias[v].dims]\n", + "df_annual_temp_bias = ds_temp_bias[sel_vars].to_dataframe()\n", + "\n", + "temp_bias_calib" + ] + }, + { + "cell_type": "markdown", + "id": "a48b2676", + "metadata": {}, + "source": [ + "Now both `melt_f` and `prcp_fac` are fixed to their default values. These are both on the lower-end to what they have been calibrated to previously when `temp_bias` is fixed. Since both of these parameters are slightly lower, the `temp_bias` is now calibrated to be higher than 0 to match the average reference mass-balance." + ] + }, + { + "cell_type": "markdown", + "id": "5b1b8ddf", + "metadata": {}, + "source": [ + "### Calibrating permutations of two parameters:\n", + "\n", + "In `mb_calibration_from_scalar_mb` in the OGGM, calibrating two parameters works as the following:\n", + "1. Set calibrate_param1 and calibrate_param2.\n", + "2. The method fixes the third parameter that has not been set.\n", + "3. The method then fixes the parameter set to calibrate_param2 and calibrates calibrate_param1.\n", + "4. If the parameter set to calibrate_param1 can find a solution with just the one parameter, the calibration will stop.\n", + "5. If the parameter set to calibrate_param1 cannot find a solution with just one parameter alone, the method will fix the parameter to either its upper or lower limit (whichever is closer to a solution) and then move onto calibrating the parameter set to calibrate_param2.\n", + "6. If a solution is found with the second calibration attempt, the calibration will stop.\n", + "7. If no solution is found with either parameter calibrated, an error will raise and the calibration will stop with no solution found." + ] + }, + { + "cell_type": "markdown", + "id": "aabbc505", + "metadata": {}, + "source": [ + "Now we will be using this technique to calibrate two parameters.\n", + "\n", + "Now calibrating `melt_f` and `prcp_fac`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15e7b69f", + "metadata": {}, + "outputs": [], + "source": [ + "# Calibrating the melt_f and prcp_fac parameters together\n", + "mf_pf_calibration = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " calibrate_param1='melt_f',\n", + " calibrate_param2='prcp_fac',\n", + " overwrite_gdir=True)\n", + "\n", + "\n", + "file_id = '_mf_pf'\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['mf_pf_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_mf_pf:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_mf_pf = ds_mf_pf.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated melt_f parameter\n", + "sel_vars = [v for v in ds_mf_pf.variables if 'month_2d' not in ds_mf_pf[v].dims]\n", + "df_annual_mf_pf = ds_mf_pf[sel_vars].to_dataframe()\n", + "\n", + "mf_pf_calibration" + ] + }, + { + "cell_type": "markdown", + "id": "d011eaae", + "metadata": {}, + "source": [ + "Notice here that we fix `temp_bias` to the default value of 0. The way that scalar calibration works, `melt_f` is calibrated first to see if a solution can be found. Therefore when the parameters are calibrated, the `prcp_fac` is firstly fixed while `melt_f` is calibrated. Since a solution can be found by just calibrating `melt_f`, as we saw at the start of this tutorial. This is why this result is the same as the calibration from just calibrating the `melt_f`." + ] + }, + { + "cell_type": "markdown", + "id": "7377e76a", + "metadata": {}, + "source": [ + "Now calibrating `prcp_fac` and `temp_bias`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6a6ef2c", + "metadata": {}, + "outputs": [], + "source": [ + "# Calibrating the prcp_fac and temp_bias parameters together\n", + "pf_tb_calibration = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " calibrate_param1='prcp_fac',\n", + " calibrate_param2='temp_bias',\n", + " overwrite_gdir=True)\n", + "\n", + "file_id = '_pf_tb'\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['pf_tb_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_pf_tb:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_pf_tb = ds_pf_tb.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated prcp_fac parameter\n", + "sel_vars = [v for v in ds_pf_tb.variables if 'month_2d' not in ds_pf_tb[v].dims]\n", + "df_annual_pf_tb = ds_pf_tb[sel_vars].to_dataframe()\n", + "\n", + "pf_tb_calibration" + ] + }, + { + "cell_type": "markdown", + "id": "dd9d286c", + "metadata": {}, + "source": [ + "Notice here that we fix `melt_f` to the default value of 5.0. To maintain the same average mass balance as before, `prcp_fac` is reduced to compensate. The parameters `melt_f` and `prcp_fac` are showing mass leaving and entering the glacier respectively." + ] + }, + { + "cell_type": "markdown", + "id": "40aaa9a3", + "metadata": {}, + "source": [ + "Now calibrating `melt_f` and `temp_bias`, we fix the `prcp_fac`. To calibrate we firstly fix the `temp_bias` and calibrate the `melt_f`, then again since a solution can be found by calibrating `melt_f` alone, the `temp_bias` does not get calibrated here:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6cae2e59", + "metadata": {}, + "outputs": [], + "source": [ + "# Calibrating the melt_f and temp_bias parameters together\n", + "mf_tb_calib = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " calibrate_param1='melt_f',\n", + " calibrate_param2='temp_bias',\n", + " overwrite_gdir=True)\n", + "\n", + "file_id = '_mf_tb'\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['mf_tb_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_mf_tb:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_mf_tb = ds_mf_tb.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated melt_f parameter\n", + "sel_vars = [v for v in ds_mf_tb.variables if 'month_2d' not in ds_mf_tb[v].dims]\n", + "df_annual_mf_tb = ds_mf_tb[sel_vars].to_dataframe()\n", + "\n", + "mf_tb_calib" + ] + }, + { + "cell_type": "markdown", + "id": "246378e3", + "metadata": {}, + "source": [ + "Now calibrating `melt_f`, `prcp_fac` and `temp_bias`, we firstly calibrate `melt_f` and fix the remaining two parameters to the default values. Since we know from previous calibrations, with the values chosen for our calibration, a solution can be found with just calibrating `melt_f` alone and therefore `prcp_fac` and `temp_bias` do not reach calibration. Therefore returning the same results as before in this tutorial where `melt_f` is calibrated first:" + ] + }, + { + "cell_type": "markdown", + "id": "d58ddcbf", + "metadata": {}, + "source": [ + "### Calibrating all three parameters:\n", + "\n", + "Now we will use `mb_calibration_from_scalar_mb` to calibrate the three parameters. This works as the following:\n", + "1. Set calibrate_param1, calibrate_param2 and calibrate_param3.\n", + "2. Fix the parameters assigned to calibrate_param2 and calibrate_param3.\n", + "3. Calibrate calibrate_param1.\n", + "4. If a solution is found, stop calibration here.\n", + "5. If no solution can be found, move onto calibrate the next parameter.\n", + "6. Repeat steps 4 and 5, therefore either reaching a calibrated solution or, if no solution can be found for any parameter, an error is raised and calibration stops with no solution.\n", + "\n", + "Here we can see that the chosen order to calibrate parameters can affect the outcomes of the calibration." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7077bbef", + "metadata": {}, + "outputs": [], + "source": [ + "# Calibrate all of the parameters\n", + "calib_all = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " calibrate_param1='melt_f',\n", + " calibrate_param2='prcp_fac',\n", + " calibrate_param3='temp_bias',\n", + " overwrite_gdir=True)\n", + "\n", + "file_id = '_calib_all'\n", + "\n", + "# Adding the mass balance model with the new calibration\n", + "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", + "fls = gdir_hef.read_pickle('inversion_flowlines')\n", + "mbdf['calib_all_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\n", + "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "# Run this again with the calibrated parameters\n", + "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_calib:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_calib = ds_calib.isel(time=slice(0, -1)).load()\n", + "\n", + "# Plot the runoff again for the calibrated melt_f parameter\n", + "sel_vars = [v for v in ds_calib.variables if 'month_2d' not in ds_calib[v].dims]\n", + "df_annual_calib = ds_calib[sel_vars].to_dataframe()\n", + "\n", + "calib_all" + ] + }, + { + "cell_type": "markdown", + "id": "b7e6bb8f", + "metadata": {}, + "source": [ + "We can see above that the calibrated parameters are now the same as in cases where we calibrate the `melt_f` first, and when we calibrate `melt_f` and `prcp_fac`. This is due to the behaviour of the scalar calibration function, which calibrates one parameter at a time in the order that we pass into the calibration function.\n", + "\n", + "This also means that the order in which we calibrate the parameters will change the outcome of our calibration, so keep this in mind when using the `mb_calibration_from_scalar_mb` function. In this tutorial we have kept to the same order of calibration in the permutation:\n", + "1. `melt_f`\n", + "2. `prcp_fac`\n", + "3. `temp_bias`\n", + "\n", + "But this order can be changed and is up to the discretion of the user to decide what they would like to calibrate initially." + ] + }, + { + "cell_type": "markdown", + "id": "31fa073d", + "metadata": {}, + "source": [ + "Now lets investigate the effect of the parameter calibration on the runoff, the aim of this tutorial!\n", + "\n", + "First we will convert the runoff values to megatonnes for readability." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0288677b", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert them to megatonnes (instead of kg)\n", + "df_runoff_melt_f = df_annual_melt_f[runoff_vars] * 1e-9\n", + "df_runoff_prcp_fac = df_annual_prcp_fac[runoff_vars] * 1e-9\n", + "df_runoff_temp_bias = df_annual_temp_bias[runoff_vars] * 1e-9\n", + "df_calib = df_annual_calib[runoff_vars] * 1e-9\n", + "df_runoff_mf_pf = df_annual_mf_pf[runoff_vars] * 1e-9\n", + "df_runoff_pf_tb = df_annual_pf_tb[runoff_vars] * 1e-9\n", + "df_runoff_mf_tb = df_annual_mf_tb[runoff_vars] * 1e-9" + ] + }, + { + "cell_type": "markdown", + "id": "1c34f50d", + "metadata": {}, + "source": [ + "And plotting runoff output values resulting from all calibration of parameters..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87e51a4a", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 3.5), sharex=True)\n", + "df_calib.sum(axis=1).plot(ax=ax, label='all params calibrated');\n", + "df_runoff_melt_f.sum(axis=1).plot(ax=ax, label='calibrated melt_f');\n", + "df_runoff.sum(axis=1).plot(ax=ax, label='default calib');\n", + "df_runoff_prcp_fac.sum(axis=1).plot(ax=ax, label='calibrated prcp_fac');\n", + "df_runoff_temp_bias.sum(axis=1).plot(ax=ax, label='calibrated temp_bias');\n", + "df_runoff_mf_pf.sum(axis=1).plot(ax=ax, label='melt_f & prcp_fac calibrated', linestyle='dashed');\n", + "df_runoff_pf_tb.sum(axis=1).plot(ax=ax, label='prcp_fac & temp_bias calibrated', linestyle='dashed');\n", + "df_runoff_mf_tb.sum(axis=1).plot(ax=ax, label='melt_f & temp_bias calibrated', linestyle='dashed');\n", + "plt.ylabel('Runoff in Megatonnes'); plt.xlabel('Years'); plt.title(f'Total annual runoff for {rgi_id}');\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "5955278f", + "metadata": {}, + "source": [ + "The impact of calibrating different permutations of parameters can be seen above, this seems to impact the runoff result quite significantly depending on the calibration choices made! It can be seen that the rough shape is similar for all calibrations, but that different permutation of calibrating parameters can have a large impact on the amount of runoff predicted by the OGGM.\n", + "\n", + "We can see that some of these runoff values are the same for certain calibrations, note that here that this is due to the style of calibration of our parameters, as we can see above the calibrations that calibrate `prcp_fac` first have a significantly lower runoff than the rest." + ] + }, + { + "cell_type": "markdown", + "id": "8bad322a", + "metadata": {}, + "source": [ + "Now plotting the mass-balance results together and comparing this to the runoff outputs:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9455a05", + "metadata": {}, + "outputs": [], + "source": [ + "# Now plotting all the mass balance results together\n", + "\n", + "# Adding the in-situ mass balance observations for comparison\n", + "mbdf['in_situ_mb'] = gdir_hef.get_ref_mb_data().loc[2000:2019]['ANNUAL_BALANCE']\n", + "\n", + "# The mass-balance from diffrent calibration runs throughout the tutorial\n", + "plt.plot(mbdf['in_situ_mb'], label='In-situ MB')\n", + "plt.plot(mbdf['melt_f_mb'], label='melt_f')\n", + "plt.plot(mbdf['prcp_fac_mb'], label='prcp_fac')\n", + "plt.plot(mbdf['temp_bias_mb'], label='temp_bias')\n", + "plt.plot(mbdf['calib_all_mb'], label='All params', linestyle='dashed')\n", + "plt.plot(mbdf['mf_pf_mb'], label='melt_f & prcp_fac', linestyle='dashed')\n", + "plt.plot(mbdf['pf_tb_mb'], label='prcp_fac & temp_bias', linestyle='dashed')\n", + "plt.plot(mbdf['mf_tb_mb'], label='melt_f & temp_bias', linestyle='dashed')\n", + "plt.legend()\n", + "plt.xlabel(\"Year\")\n", + "plt.ylabel(\"Mass Balance\");\n", + "plt.title(f'Mass Balance time series for {rgi_id}');" + ] + }, + { + "cell_type": "markdown", + "id": "e4aa15ef", + "metadata": {}, + "source": [ + "In this graph we can now see that our different calibrations also have an impact on the mass-balance\n", + "\n", + "Here, we can see that there are a few overlaps again for different permutations of calibrations, if we look closely we can see that these are the same as the above for runoff.\n", + "\n", + "However, from initial inspection we can see that for the runoff there is quite a significant difference based off of choices of calibrations. Here we can see that the differences between the mass-balance don't seem quite as dramatic.\n", + "\n", + "From the runoff graph, the calibration of the `prcp_fac` appears to be impactful on the runoff. We can now explore a bit further, and investigate the relationship between `prcp_fac` and `melt_on_glacier`, which is a component of the `runoff`. We will do this by fixing the `prcp_fac` to a range of values and calibrating the `melt_f` to obtain a range of calibrated values. This will allow us to verify the sensitivity of these parameters." + ] + }, + { + "cell_type": "markdown", + "id": "eef65f7e", + "metadata": {}, + "source": [ + "## Sensitivity Analysis of runoff parameters in calibration\n", + "\n", + "**Sensitivity Analysis** is a technique used to determine how a model output is affected by the changes in parameters. We will perform a simple sensitvity analysis in this tutorial to investigate the affects of parameter calibration on some of the runoff components.\n", + "\n", + "We will now investigate the sensitivity of the runoff parameters to the calibration of the mass-balance parameters below. This will allow us to understand the relationship between parameters further, for example if one parameter is changed, how much does this affect the other parameter?" + ] + }, + { + "cell_type": "markdown", + "id": "af1a9529", + "metadata": {}, + "source": [ + "**First let's start with some definitions!**\n", + "\n", + "We will be investigating the melt contribution, which is:\n", + "$ \\frac{\\text{melt on glacier}}{\\text{runoff`}}$" + ] + }, + { + "cell_type": "markdown", + "id": "003110aa", + "metadata": {}, + "source": [ + "In this tutorial we will investigate the sensitivities of the `melt_on_glacier` and melt contribution to the `prcp_fac`.\n", + "\n", + "We have chosen these parameters as `melt_on_glacier` is a single variable, defined by a linear relationship from the `prcp_fac`. Here we will just be affecting the `melt_on_glacier` directly.\n", + "\n", + "The melt contribution variable is a bit more complicated, as this will be investigating the relationship between the run off which includes `melt_on_glacier` variable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "834456bf", + "metadata": {}, + "outputs": [], + "source": [ + "# Varying prcp_fac between a range of values with a step of 1.0\n", + "pd_prcp_sens = pd.DataFrame(index=np.arange(0.1, 5.0, 0.3))\n", + "file_id = '_sens'\n", + "\n", + "# Calibrate the melt factor for each precipitation factor\n", + "spec_mb_melt_f_sens_dict = {}\n", + "for pf in pd_prcp_sens.index:\n", + " calib_param = mb_calibration_from_scalar_mb(gdir_hef,\n", + " ref_mb = ref_mb,\n", + " prcp_fac = pf,\n", + " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", + " overwrite_gdir=True)\n", + "\n", + " # Fill the dataframe with the calibrated parameters\n", + " pd_prcp_sens.loc[pf, 'melt_f'] = calib_param['melt_f']\n", + "\n", + " # We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + " # Run this again with the calibrated parameters\n", + " tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "\n", + " with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", + "\n", + " # Plot the runoff again for the calibrated melt_f parameter\n", + " sel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + " df_annual_sens = ds_sens[sel_vars].to_dataframe()\n", + "\n", + " pd_prcp_sens.loc[pf, 'melt_off_glacier'] = df_annual_sens['melt_off_glacier'].mean()\n", + " pd_prcp_sens.loc[pf, 'melt_on_glacier'] = df_annual_sens['melt_on_glacier'].mean()\n", + " pd_prcp_sens.loc[pf, 'liq_prcp_off_glacier'] = df_annual_sens['liq_prcp_off_glacier'].mean()\n", + " pd_prcp_sens.loc[pf, 'liq_prcp_on_glacier'] = df_annual_sens['liq_prcp_on_glacier'].mean()\n", + " pd_prcp_sens.loc[pf, 'runoff'] = pd_prcp_sens.loc[pf, 'melt_off_glacier'] + pd_prcp_sens.loc[pf, 'melt_on_glacier'] + pd_prcp_sens.loc[pf, 'liq_prcp_off_glacier'] + pd_prcp_sens.loc[pf, 'liq_prcp_on_glacier']\n", + "\n", + "colors_melt_f = plt.get_cmap('viridis').colors[10::10]\n", + "plt.figure()\n", + "for j, pf in enumerate(pd_prcp_sens.index):\n", + " plt.plot(pf, pd_prcp_sens.loc[pf, 'melt_on_glacier'], 'o', color=colors_melt_f[j])\n", + " plt.xlabel('Precipitation factor')\n", + " plt.ylabel('Mean annual melt on glacier (kg m-2 yr-1)')\n", + " plt.title('Sensitivity of Melt On Glacier to Precipitation Factor')\n", + "\n", + "# Here we are varying the precipitation factor and calibrating the melt factor for each value\n", + "# Then plotting the mean annual melt_on_glacier against the precipitation factor" + ] + }, + { + "cell_type": "markdown", + "id": "00bed96c", + "metadata": {}, + "source": [ + "In the above plot, a strong linear positive correlation is shown between the precipitation factor and the mean annual melt on the glacier. \n", + "\n", + "From the calibration, it can be seen that as the precipitation factor is changed, the melt factor parameter is calibrated accordingly and increases with the precipitation factor to maintain the average annual mass-balance. Therefore as the `prcp_fac` increases, the `melt_on_glacier` increases linearly in a 1:1 ratio.\n", + "\n", + "In the OGGM, the `melt_on_glacier` is derived from both the `melt_f` parameter and `prcp_fac`. Therefore this sensitivity makes sense." + ] + }, + { + "cell_type": "markdown", + "id": "28875fe5", + "metadata": {}, + "source": [ + "Now investigating the sensitivity of the Glacier Melt Contribution to the Precipitation Factor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ad8fa07", + "metadata": {}, + "outputs": [], + "source": [ + "colors_melt_f = plt.get_cmap('viridis').colors[10::10]\n", + "plt.figure()\n", + "for j, pf in enumerate(pd_prcp_sens.index):\n", + " plt.plot(pf, pd_prcp_sens.loc[pf, 'melt_on_glacier']/pd_prcp_sens.loc[pf, 'runoff'], 'o', color=colors_melt_f[j])\n", + " plt.xlabel('Precipitation factor')\n", + " plt.ylabel('Mean Glacial Melt Contribution')\n", + " plt.title('Sensitivity of Glacial Melt Contribution to the Precipitation Factor')\n", + "\n", + "# Here we are plotting how the glacial melt contribution (melt_on_glacier/runoff) varies with the changing precipitation factor\n", + "# This shows how the relative contribution of glacier melt to total runoff changes as we adjust precipitation" + ] + }, + { + "cell_type": "markdown", + "id": "75444665", + "metadata": {}, + "source": [ + "It can be seen that there is a relatively strong negative correlation, between the precipitation factor and the mean glacial melt contribution. And therefore it can be seen that the mean glacial melt contribution is sensitive to the change in precipitation factor. As we increase the `prcp_fac`, we know that the `melt_f` increases along with this. Considering the melt contribution variable contains a numerator of `melt_on_glacier`, which has been shown to increase as the `prcp_fac` parameter increases.\n", + "\n", + "The denominator is the runoff, which contains 4 summed variables (including `melt_on_glacier`) that are affected by the `prcp_fac`, all of which increase as the `prcp_fac` increases, as can be seen within the OGGM model implementation. The runoff therefore increases at least as much as the `melt_on_glacier`.\n", + "\n", + "Therefore the increase of the denominator overpowers the increase of the numerator, and the relationship shown above reflects this reaction to the increase in `prcp_fac`." + ] + }, + { + "cell_type": "markdown", + "id": "cf62afed", + "metadata": {}, + "source": [ + "## Exercise:\n", + "\n", + "Try this yourself for sensitivity with the precipitation factor when varying the melt factor as above:\n", + "1. The melt_off_glacier parameter.\n", + "2. The total liquid precipitation parameters (sum of on and off).\n", + "3. Any other interesting combinations you can discover!\n", + "\n", + "Does anything you discover surprise you?\n", + "Are there particularly sensitive parameters compared to others?\n", + "\n", + "**Stretch:** Can you experiment different combinations of parameters for sensitvity, e.g. such as varying the temp_bias parameters instead? " + ] + }, + { + "cell_type": "markdown", + "id": "8e0eae95", + "metadata": {}, + "source": [ + "## WGMS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53b9a503", + "metadata": {}, + "outputs": [], + "source": [ + "# # We start by fetching the observations\n", + "# mbdf = gdir_hef.get_ref_mb_data()[['ANNUAL_BALANCE']]\n", + "# mbdf.columns = ['in_situ_mb']\n", + "# mbdf = mbdf.loc[2000:2019]\n", + "# mbdf.plot();\n", + "# mbdf.mean()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "155dbd5a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/welcome.ipynb b/notebooks/welcome.ipynb index 4a2399de..61a349a1 100644 --- a/notebooks/welcome.ipynb +++ b/notebooks/welcome.ipynb @@ -76,7 +76,12 @@ "source": [ "## Hydrological output\n", "\n", - "- [Hydrological mass-balance output](tutorials/hydrological_output.ipynb)" + "- [Hydrological mass-balance output](tutorials/hydrological_output.ipynb)\n", + "- [Sensitivity of glacier runoff to the mass balance parameters](tutorials/runoff_sensitivity.ipynb)\n", + "\n", + "You might find the following notebooks in OGGM-Edu interesting as well!\n", + "- [Glaciers as water resources: part 1 (idealized climate)](https://oggm.org/oggm-edu-notebooks/oggm-edu/glacier_water_resources.html)\n", + "- [Glaciers as water resources: part 2 (projections)](https://oggm.org/oggm-edu-notebooks/oggm-edu/glacier_water_resources_projections.html)" ] }, { From 2b008345501e538a99b2d96231a8d96f2393680d Mon Sep 17 00:00:00 2001 From: Chloe Hancock Date: Tue, 7 Apr 2026 09:58:00 +0100 Subject: [PATCH 2/5] Removed uneeded code in sensitivity notebook --- notebooks/tutorials/runoff_sensitivity.ipynb | 73 +------------------- 1 file changed, 3 insertions(+), 70 deletions(-) diff --git a/notebooks/tutorials/runoff_sensitivity.ipynb b/notebooks/tutorials/runoff_sensitivity.ipynb index 0184422f..8ccae4a8 100644 --- a/notebooks/tutorials/runoff_sensitivity.ipynb +++ b/notebooks/tutorials/runoff_sensitivity.ipynb @@ -187,75 +187,6 @@ "### **Calibration:** Scalar Mass Balance Calibration" ] }, - { - "cell_type": "markdown", - "id": "d5f9d1ac", - "metadata": {}, - "source": [ - "Firstly comparing the mass balances below for in-situ observations and the default calibrated parameters. To understand the current behaviour of the mass balance output.\n" - ] - }, - { - "cell_type": "markdown", - "id": "fe1bd5be", - "metadata": {}, - "source": [ - "Now let's experiment by calibrating the parameters: these parameters are the melt factor (`melt_f`), the precipitation factor (`prcp_fac`) and the temperature bias (`temp_bias`). We will calibrate each parameter and permutations of the calibrations:\n", - "1. Just calibrating `melt_f`\n", - "2. Just calibrating `prcp_fac`\n", - "3. Just calibrating `temp_bias`\n", - "4. Calibrating `melt_f` and `prcp_fac`\n", - "5. Calibrating `prcp_fac` and `temp_bias`\n", - "6. Calibrating `melt_f` and `temp_bias`\n", - "7. Calibrating all parameters; `melt_f`, `prcp_fac` and `temp_bias`\n", - "\n", - "The next calibration section will be relatively repetative, but take note of the parameters being calibrated and the inputs being used in the calibration process.\n", - "\n", - "We will start the calibration by just calibrating the model the `melt_f` paramater." - ] - }, - { - "cell_type": "markdown", - "id": "33360649", - "metadata": {}, - "source": [ - "
\n", - " \n", - " Note: In the OGGM Scalar Calibration method, even when all parameters are nominated for calibration, they may not reach the calibration stage. The OGGM works by calibrating one parameter at a time, only moving onto the next when the previous parameter calibration cannot reach the desired target. In the tutorial below, we will see the result of this in action with usually only the first parameter being calibrated, even if all parameters are nominated for calibration.\n", - " \n", - "
" - ] - }, - { - "cell_type": "markdown", - "id": "f47b6995", - "metadata": {}, - "source": [ - "### Single Parameter Calibration" - ] - }, - { - "cell_type": "markdown", - "id": "d28b4228", - "metadata": {}, - "source": [ - "Fetch the reference mass balance we want to calibrate to from the geodetic mass-balance data, this is the data from Hugonnet et al 2021:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0e675125", - "metadata": {}, - "outputs": [], - "source": [ - "ref_mb_df = utils.get_geodetic_mb_dataframe().loc[gdir_hef.rgi_id]\n", - "ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == cfg.PARAMS['geodetic_mb_period']].iloc[0]\n", - "# dmdtda: in meters water-equivalent per year -> we convert to kg m-2 yr-1\n", - "ref_mb = ref_mb_df['dmdtda'] * 1000\n", - "ref_mb" - ] - }, { "cell_type": "markdown", "id": "3ec4b550", @@ -271,6 +202,8 @@ "metadata": {}, "outputs": [], "source": [ + "ref_mb = -1100\n", + "\n", "# Just calibrate the melt_f parameter\n", "calib_param_melt_f = mb_calibration_from_scalar_mb(gdir_hef,\n", " ref_mb = ref_mb,\n", @@ -1034,7 +967,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.0" + "version": "3.11.14" } }, "nbformat": 4, From 9e030aacf1ead374ae79aa4279e4ce49800871e1 Mon Sep 17 00:00:00 2001 From: Chloe Hancock Date: Thu, 9 Apr 2026 17:15:40 +0100 Subject: [PATCH 3/5] Updated runoff sensitivity tutorial, added new plots, tightened phrasing --- notebooks/tutorials/runoff_sensitivity.ipynb | 974 +++++++------------ 1 file changed, 337 insertions(+), 637 deletions(-) diff --git a/notebooks/tutorials/runoff_sensitivity.ipynb b/notebooks/tutorials/runoff_sensitivity.ipynb index 8ccae4a8..edb788e5 100644 --- a/notebooks/tutorials/runoff_sensitivity.ipynb +++ b/notebooks/tutorials/runoff_sensitivity.ipynb @@ -1,5 +1,33 @@ { "cells": [ + { + "cell_type": "markdown", + "id": "349a2484", + "metadata": {}, + "source": [ + "# Runoff Sensitivity to Mass Balance Parameters" + ] + }, + { + "cell_type": "markdown", + "id": "0f9dec41", + "metadata": {}, + "source": [ + "Goals of this notebook: To explore the sensitivity of the runoff and OGGM's hydrological components to the mass balance parameters." + ] + }, + { + "cell_type": "markdown", + "id": "543e58a1", + "metadata": {}, + "source": [ + "Our other hydrology tutorials introduce the idea of hydrology in the OGGM. It is recommended to consult the previous tutorials to gain an understanding first as to how hydrology works in the OGGM. \n", + "\n", + "[Glaciers as water resources: Part 1](https://edu-notebooks.oggm.org/oggm-edu/glacier_water_resources.html)\n", + "\n", + "[Glaciers as water resources: Part 2](https://edu-notebooks.oggm.org/oggm-edu/glacier_water_resources_projections.html)" + ] + }, { "cell_type": "markdown", "id": "d1abcbfd", @@ -7,9 +35,7 @@ "source": [ "## Set Up\n", "\n", - "First install required packages to run this tutorial. \n", - "\n", - "We will be using level 4 data to allow for the dynamical spinup and the calibration capabilities." + "First install required packages to run this tutorial and initialize our glacier directories!" ] }, { @@ -23,10 +49,11 @@ "import pandas as pd\n", "import numpy as np\n", "import xarray as xr\n", + "import matplotlib\n", "\n", - "from oggm import cfg, utils, workflow, tasks\n", + "from oggm import cfg, utils, workflow, tasks, DEFAULT_BASE_URL\n", "from oggm.core import massbalance\n", - "from oggm.core.massbalance import mb_calibration_from_scalar_mb" + "from oggm.core.massbalance import MultipleFlowlineMassBalance" ] }, { @@ -60,764 +87,419 @@ "rgi_id = 'RGI60-11.00897'\n", "\n", "# We pick the elevation-bands glaciers because they run a bit faster - but they create more step changes in the area outputs\n", - "from oggm import DEFAULT_BASE_URL\n", - "gdir_hef = workflow.init_glacier_directories([rgi_id], from_prepro_level=4, prepro_border=80, prepro_base_url=DEFAULT_BASE_URL)[0]" + "gdir_hef = workflow.init_glacier_directories([rgi_id], from_prepro_level=4, prepro_border=160, prepro_base_url=DEFAULT_BASE_URL)[0]" ] }, { "cell_type": "markdown", - "id": "78de2a38", - "metadata": {}, - "source": [ - "## Generating the Hydrological and Glaciological outputs" - ] - }, - { - "cell_type": "markdown", - "id": "041261fe", - "metadata": {}, - "source": [ - "The pre-processed directories don't automatically have hydrological outputs, so lets run the `run_with_hydro` task below to calculate these! This requires a dynamical spinup, let's choose the time period 2000-2020." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "36a47e5d", + "id": "eef65f7e", "metadata": {}, - "outputs": [], "source": [ - "# file identifier where the model output is saved\n", - "file_id = '_default'\n", + "# An Introduction to Sensitivity Analysis\n", "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # which climate scenario? See following notebook for other examples\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "**Sensitivity Analysis** investigates how the variation in the output of a numerical model can be attributed to variations of its input factors [1](https://www.sciencedirect.com/science/article/pii/S1364815216300287).\n", "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix='_default')) as ds:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds = ds.isel(time=slice(0, -1)).load()" - ] - }, - { - "cell_type": "markdown", - "id": "ac645b31", - "metadata": {}, - "source": [ - "Now that we have the hydrological and glaciological outputs for our desired time period, lets calculate the total runoff for the glacier. There are four key variables that form the total runoff when summed together, and the process to calculate the runoff can be seen below.\n", + "In this tutorial, we perform a simple, exploratory one-at-a-time sensitivity analysis to investigate the effects of the mass balance parameters on the glaciohydrological model outputs. We will focus on the mass balance parameters: the melt factor, temperature bias and precipitation factor.\n", "\n", - "The total runoff from the glacier is:\n", - "`runoff` = `melt_off_glacier` + `melt_on_glacier` + `liq_prcp_off_glacier` + `liq_prcp_on_glacier`\n" + "This will allow us to understand the relationship between our parameters and model output, for example if one parameter is changed, how much does this affect our model output?" ] }, { "cell_type": "markdown", - "id": "1d20120d", - "metadata": {}, - "source": [ - "The figure below shows the runoff for the default calibrated parameters." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d27a8a7d", + "id": "bb462087", "metadata": {}, - "outputs": [], "source": [ - "sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]\n", - "df_annual = ds[sel_vars].to_dataframe()\n", - "\n", - "# These summed variabels give the total runoff from the glacier\n", - "runoff_vars = ['melt_off_glacier', 'melt_on_glacier','liq_prcp_off_glacier', 'liq_prcp_on_glacier']\n", - "\n", - "# Convert them to megatonnes (instead of kg)\n", - "df_runoff = df_annual[runoff_vars] * 1e-9\n", - "fig, ax = plt.subplots(figsize=(10, 3.5), sharex=True)\n", - "df_runoff.sum(axis=1).plot(ax=ax);\n", - "plt.ylabel('Runoff'); plt.xlabel('Years'); plt.title(f'Total annual runoff for {rgi_id}');\n" + "# Simple Sensitivity Analysis Experiment" ] }, { "cell_type": "markdown", - "id": "92501cfa", + "id": "050f0271", "metadata": {}, "source": [ - "This was done with the following parameter set:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0a217ac7", - "metadata": {}, - "outputs": [], - "source": [ - "gdir_hef.read_json('mb_calib')" + "### Temperature Bias" ] }, { "cell_type": "markdown", - "id": "33a812b6", + "id": "a9c0ce89", "metadata": {}, "source": [ - "## Parameter Calibration" - ] - }, - { - "cell_type": "markdown", - "id": "0d96ff15", - "metadata": {}, - "source": [ - "The mass-balance parameters at Level 4 are pre-calibrated automatically by the OGGM and stored in the dataframe. However there are flexible mass-balance calibration schemes in the OGGM.\n", - "\n", - "In this tutorial, our aim is to calibrate the mass-balance parameters and investigate the impact that this calibration has on the run-off. In the OGGM we can currently (at v1.6.3) calibrate the `melt_f`, `prcp_fac` and `temp_bias` parameters. We will calibrate these parameters using the **scalar mass balance** method, and investigate the relationship with the runoff components and the total runoff output." - ] - }, - { - "cell_type": "markdown", - "id": "ac6c6eaa", - "metadata": {}, - "source": [ - "### **Calibration:** Scalar Mass Balance Calibration" - ] - }, - { - "cell_type": "markdown", - "id": "3ec4b550", - "metadata": {}, - "source": [ - "Now we will calibrate with the scalar mass balance and use its default calibration setting. This is just calibrating the `melt_f` parameter:" + "We will begin by investigating the sensitvity of the runoff to one parameter at a time, we will start with **temperature bias**. Below, we will vary only the temperature bias, and fix the melt factor and the precipitation factor." ] }, { "cell_type": "code", "execution_count": null, - "id": "5b5dab26", + "id": "2d1cf29c", "metadata": {}, "outputs": [], "source": [ - "ref_mb = -1100\n", + "temp_bias_dict = {}\n", + "file_id = '_sens'\n", "\n", - "# Just calibrate the melt_f parameter\n", - "calib_param_melt_f = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " overwrite_gdir=True)\n", + "for temp_bias in np.arange(-5,5.0,0.5): # We are varying the temperature bias\n", "\n", - "# Creating a mass balance data frame to store the results\n", - "mbdf= pd.DataFrame(index = np.arange(2000,2020,1))\n", + "\t# For each temperature bias, we create a mass balance model with the same melt factor and precipitation factor, but with the new temperature bias\n", + "\tmb_model = MultipleFlowlineMassBalance(\n", + " \tgdir_hef,\n", + " \tmb_model_class=massbalance.MonthlyTIModel,\n", + " \ttemp_bias=float(temp_bias), # Vary the temperature bias\n", + "\t\tmelt_f=5.0, # Fix melt factor to 5.0 for all runs, to see the effect of the temperature bias alone\n", + "\t\tprcp_fac=2.5, # Fix precipitation factor to 2.5 for all runs, to see the effect of the temperature bias alone\n", + " \tcheck_calib_params=False,\n", + " \t) \n", + "\t\n", + "\t# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "\t# Run this with the our 2 fixed mass balance parameters and our varying temperature bias\n", + "\ttasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + "\t\t\t\t \t\trun_task=tasks.run_from_climate_data, # Run from climate data\n", + "\t\t\t\t \t\tmb_model=mb_model, # use the mass balance model with the new temp_bias\n", + "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", + "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", + "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", + "\t\n", + "\t# Now read the hydrological outputs for this run\n", + "\twith xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + "\t\t# The last step of hydrological output is NaN (we can't compute it for this year)\n", + "\t\tds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['melt_f_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", "\n", - "calib_param_melt_f" + "\t# Calcuating the runoff for the temperature bias value and storing it in a dictionary to plot\n", + "\ttemp_bias_dict[temp_bias] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier'])*1e-9 # Convert from Kilograms to Megatonnes per year (Mt/yr) for easier plotting" ] }, { "cell_type": "markdown", - "id": "b2f0fcc5", + "id": "96c586d6", "metadata": {}, "source": [ - "We can see the parameters above, we will use the run_with_hydro task to gather the hydrological outputs. We will now calculate the runoff from the calabrated parameters." + "Now, for each temperature bias value, let's plot the mean runoff against the temperature bias, and plot the runoff time series to undestand how this varies across our range of temperature bias values." ] }, { "cell_type": "code", "execution_count": null, - "id": "97847d11", + "id": "801bc112", "metadata": {}, "outputs": [], "source": [ - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "file_id = '_melt_f'\n", + "# Now let's get a nice colormap centered at temp_bias=0\n", + "norm = matplotlib.colors.Normalize(vmin=-5, vmax=5.01)\n", + "colors_temp_bias = plt.get_cmap('coolwarm')\n", "\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", + "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", + "for j, temp_bias in enumerate(temp_bias_dict.keys()):\n", + " axs[0].plot(temp_bias, temp_bias_dict[temp_bias].mean(), 'o',\n", + " color=colors_temp_bias(norm(temp_bias))) \n", + "axs[0].set_ylabel('Mean Runoff (Mt/yr)')\n", + "axs[0].set_xlabel('temp_bias (°C)');\n", + "axs[0].set_title('Mean Runoff vs Temperature Bias')\n", "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_melt_f:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_melt_f = ds_melt_f.isel(time=slice(0, -1)).load()\n", "\n", - "# Plot the runoff again for the calibrated melt_f parameter\n", - "sel_vars = [v for v in ds_melt_f.variables if 'month_2d' not in ds_melt_f[v].dims]\n", - "df_annual_melt_f = ds_melt_f[sel_vars].to_dataframe()\n" + "for temp_bias in temp_bias_dict.keys():\n", + " axs[1].plot(np.arange(2000,2020,1),\n", + " temp_bias_dict[temp_bias].values, '-', \n", + " color=colors_temp_bias(norm(temp_bias)),\n", + " label=temp_bias)\n", + "axs[1].set_ylabel('Runoff (Mt/yr)')\n", + "axs[1].set_xlabel('Year')\n", + "axs[1].legend(title='temp_bias:', bbox_to_anchor=(1,1))\n", + "axs[1].set_xticks(np.arange(2000,2020,2));\n", + "axs[1].set_title('Runoff Time Series for Different Temperature Biases')\n", + "\n", + "fig.suptitle(\"Runoff and Temperature Bias - {}\".format(gdir_hef.rgi_id));\n", + "plt.tight_layout()" ] }, { "cell_type": "markdown", - "id": "2f0173ff", + "id": "e6917d8d", "metadata": {}, "source": [ - "Now calibrating the precipitation factor variable, `prcp_fac`. Here we fix the `melt_f` and the `temp_bias` and then calibrate the precipitation factor accordingly, this will ensure that the reference average mass balance will be matched." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f4182717", - "metadata": {}, - "outputs": [], - "source": [ - "# Just calibrate the prcp_fac parameter\n", - "calib_param_prcp_fac = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " calibrate_param1='prcp_fac',\n", - " overwrite_gdir=True)\n", - "\n", - "file_id = '_prcp_fac'\n", - "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['prcp_fac_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", - "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", - "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_prcp_fac:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_prcp_fac = ds_prcp_fac.isel(time=slice(0, -1)).load()\n", - "\n", - "# Plot the runoff again for the calibrated melt_f parameter\n", - "sel_vars = [v for v in ds_prcp_fac.variables if 'month_2d' not in ds_prcp_fac[v].dims]\n", - "df_annual_prcp_fac = ds_prcp_fac[sel_vars].to_dataframe()\n", - "\n", - "calib_param_prcp_fac" + "We can see that the runoff appears to be sensitive to the temperature bias! In the graph on the left, we can see the mean runoff increases as we increase the temperature bias, and on the graph on the right, we can see how this affects the runoff annually." ] }, { "cell_type": "markdown", - "id": "86551587", + "id": "e2461a95", "metadata": {}, "source": [ - "Note here that now the `melt_f` is lower than before, and the `prcp_fac` is now lowered in calibration so the average reference mass-balance can be matched." + "### Precipitation Factor" ] }, { "cell_type": "markdown", - "id": "30085bd7", + "id": "8ad1e197", "metadata": {}, "source": [ - "Now calibrating the temperature bias parameter, `temp_bias`. Here we now fix the `melt_f` and `prcp_fac` to the default values, and calibrate the `temp_bias` to match the average reference mass-balance:" + "Now lets explore what happens when we very the **precipitation factor**, and fix our other two mass balance parameters." ] }, { "cell_type": "code", "execution_count": null, - "id": "0279da19", + "id": "b659f889", "metadata": {}, "outputs": [], "source": [ - "# Just calibrate the temp_bias parameters\n", - "temp_bias_calib = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " calibrate_param1='temp_bias',\n", - " overwrite_gdir=True)\n", - "\n", - "file_id = '_temp_bias'\n", - "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['temp_bias_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", - "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", - "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_temp_bias:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_temp_bias = ds_temp_bias.isel(time=slice(0, -1)).load()\n", + "prcp_fac_dict = {}\n", "\n", - "# Plot the runoff again for the calibrated melt_f parameter\n", - "sel_vars = [v for v in ds_temp_bias.variables if 'month_2d' not in ds_temp_bias[v].dims]\n", - "df_annual_temp_bias = ds_temp_bias[sel_vars].to_dataframe()\n", - "\n", - "temp_bias_calib" - ] - }, - { - "cell_type": "markdown", - "id": "a48b2676", - "metadata": {}, - "source": [ - "Now both `melt_f` and `prcp_fac` are fixed to their default values. These are both on the lower-end to what they have been calibrated to previously when `temp_bias` is fixed. Since both of these parameters are slightly lower, the `temp_bias` is now calibrated to be higher than 0 to match the average reference mass-balance." - ] - }, - { - "cell_type": "markdown", - "id": "5b1b8ddf", - "metadata": {}, - "source": [ - "### Calibrating permutations of two parameters:\n", + "for prcp_fac in np.arange(0.1,10,0.5): # Now we are varying the precipitation factor\n", + "\t\n", + "\t# For each precipitation factor, we create a mass balance model with the same melt factor and temperature bias, but with the new precipitation factor\n", + "\tmb_model = MultipleFlowlineMassBalance(\n", + " \tgdir_hef,\n", + " \tmb_model_class=massbalance.MonthlyTIModel,\n", + " \tprcp_fac=float(prcp_fac),\n", + "\t\tmelt_f=5.0, # Fix melt factor\n", + "\t\ttemp_bias=0, # Fix the temperature bias to 0\n", + " \tcheck_calib_params=False,\n", + " \t) \n", + "\t\n", + "\t# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", + "\t# Run this with the our 2 fixed mass balance parameters and our varying precipitation factor\n", + "\ttasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + "\t\t\t\t \t\trun_task=tasks.run_from_climate_data,\n", + "\t\t\t\t \t\tmb_model=mb_model, # use the mass balance model with the new prcp_fac\n", + "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", + "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", + "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", + "\t\n", + "\t# Reading the glaciological and hydrological outputs for this run again\n", + "\twith xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + "\t\t# The last step of hydrological output is NaN (we can't compute it for this year)\n", + "\t\tds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", + "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", "\n", - "In `mb_calibration_from_scalar_mb` in the OGGM, calibrating two parameters works as the following:\n", - "1. Set calibrate_param1 and calibrate_param2.\n", - "2. The method fixes the third parameter that has not been set.\n", - "3. The method then fixes the parameter set to calibrate_param2 and calibrates calibrate_param1.\n", - "4. If the parameter set to calibrate_param1 can find a solution with just the one parameter, the calibration will stop.\n", - "5. If the parameter set to calibrate_param1 cannot find a solution with just one parameter alone, the method will fix the parameter to either its upper or lower limit (whichever is closer to a solution) and then move onto calibrating the parameter set to calibrate_param2.\n", - "6. If a solution is found with the second calibration attempt, the calibration will stop.\n", - "7. If no solution is found with either parameter calibrated, an error will raise and the calibration will stop with no solution found." + "\t# Calculating the runoff\n", + "\tprcp_fac_dict[prcp_fac] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier']) * 1e-9" ] }, { "cell_type": "markdown", - "id": "aabbc505", + "id": "ae390ad3", "metadata": {}, "source": [ - "Now we will be using this technique to calibrate two parameters.\n", - "\n", - "Now calibrating `melt_f` and `prcp_fac`:" + "Now let's plot to see our results." ] }, { "cell_type": "code", "execution_count": null, - "id": "15e7b69f", + "id": "9ff119cd", "metadata": {}, "outputs": [], "source": [ - "# Calibrating the melt_f and prcp_fac parameters together\n", - "mf_pf_calibration = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " calibrate_param1='melt_f',\n", - " calibrate_param2='prcp_fac',\n", - " overwrite_gdir=True)\n", - "\n", + "# Now we can centre the colormap around prcp_fac\n", + "norm = matplotlib.colors.Normalize(vmin=0.1, vmax=10)\n", + "colors_prcp_fac = plt.get_cmap('coolwarm')\n", "\n", - "file_id = '_mf_pf'\n", + "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", + "for j, prcp_fac in enumerate(prcp_fac_dict.keys()):\n", + " axs[0].plot(prcp_fac, prcp_fac_dict[prcp_fac].mean(), 'o',\n", + " color=colors_prcp_fac(norm(prcp_fac))) \n", + "axs[0].set_ylabel('Mean Runoff (Mt/yr)')\n", + "axs[0].set_xlabel('Precipitation factor');\n", + "axs[0].set_title('Mean Runoff vs Precipitation Factor')\n", "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['mf_pf_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "for prcp_fac in prcp_fac_dict.keys():\n", + " axs[1].plot(np.arange(2000,2020,1),\n", + " prcp_fac_dict[prcp_fac].values, '-', \n", + " color=colors_prcp_fac(norm(prcp_fac)),\n", + " label=prcp_fac)\n", + "axs[1].set_ylabel('Runoff (Mt/yr)')\n", + "axs[1].set_xlabel('Year')\n", + "axs[1].legend(title='Precipitation factor:', bbox_to_anchor=(1,1))\n", + "axs[1].set_xticks(np.arange(2000,2020,2));\n", + "axs[1].set_title('Runoff Time Series for Different Precipitation Factors')\n", "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", - "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_mf_pf:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_mf_pf = ds_mf_pf.isel(time=slice(0, -1)).load()\n", - "\n", - "# Plot the runoff again for the calibrated melt_f parameter\n", - "sel_vars = [v for v in ds_mf_pf.variables if 'month_2d' not in ds_mf_pf[v].dims]\n", - "df_annual_mf_pf = ds_mf_pf[sel_vars].to_dataframe()\n", - "\n", - "mf_pf_calibration" - ] - }, - { - "cell_type": "markdown", - "id": "d011eaae", - "metadata": {}, - "source": [ - "Notice here that we fix `temp_bias` to the default value of 0. The way that scalar calibration works, `melt_f` is calibrated first to see if a solution can be found. Therefore when the parameters are calibrated, the `prcp_fac` is firstly fixed while `melt_f` is calibrated. Since a solution can be found by just calibrating `melt_f`, as we saw at the start of this tutorial. This is why this result is the same as the calibration from just calibrating the `melt_f`." + "fig.suptitle(\"Runoff and Precipitation Factor - {}\".format(gdir_hef.rgi_id));\n", + "plt.tight_layout()" ] }, { "cell_type": "markdown", - "id": "7377e76a", - "metadata": {}, - "source": [ - "Now calibrating `prcp_fac` and `temp_bias`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f6a6ef2c", + "id": "44c2ecc9", "metadata": {}, - "outputs": [], "source": [ - "# Calibrating the prcp_fac and temp_bias parameters together\n", - "pf_tb_calibration = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " calibrate_param1='prcp_fac',\n", - " calibrate_param2='temp_bias',\n", - " overwrite_gdir=True)\n", - "\n", - "file_id = '_pf_tb'\n", - "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['pf_tb_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", - "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", - "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_pf_tb:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_pf_tb = ds_pf_tb.isel(time=slice(0, -1)).load()\n", - "\n", - "# Plot the runoff again for the calibrated prcp_fac parameter\n", - "sel_vars = [v for v in ds_pf_tb.variables if 'month_2d' not in ds_pf_tb[v].dims]\n", - "df_annual_pf_tb = ds_pf_tb[sel_vars].to_dataframe()\n", - "\n", - "pf_tb_calibration" + "We can see that again, the runoff appears sensitive to the precipitation factor! And that as the precipitation factor increases, as does the runoff." ] }, { "cell_type": "markdown", - "id": "dd9d286c", + "id": "0958f452", "metadata": {}, "source": [ - "Notice here that we fix `melt_f` to the default value of 5.0. To maintain the same average mass balance as before, `prcp_fac` is reduced to compensate. The parameters `melt_f` and `prcp_fac` are showing mass leaving and entering the glacier respectively." + "### Melt Factor" ] }, { "cell_type": "markdown", - "id": "40aaa9a3", + "id": "2764c928", "metadata": {}, "source": [ - "Now calibrating `melt_f` and `temp_bias`, we fix the `prcp_fac`. To calibrate we firstly fix the `temp_bias` and calibrate the `melt_f`, then again since a solution can be found by calibrating `melt_f` alone, the `temp_bias` does not get calibrated here:" + "Finally, let's see what happens when we alter the melt factor and fix the remaining 2 mass balance parameters!" ] }, { "cell_type": "code", "execution_count": null, - "id": "6cae2e59", + "id": "eb55224d", "metadata": {}, "outputs": [], "source": [ - "# Calibrating the melt_f and temp_bias parameters together\n", - "mf_tb_calib = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " calibrate_param1='melt_f',\n", - " calibrate_param2='temp_bias',\n", - " overwrite_gdir=True)\n", + "melt_f_dict = {}\n", "\n", - "file_id = '_mf_tb'\n", + "for melt_f in np.arange(1.5,17,1.0): # Now we are varying the melt factor\n", + "\t\n", + "\t# For each melt factor, we create a mass balance model with the same precipitation factor and temperature bias, but now with the new melt factor\n", + "\tmb_model = MultipleFlowlineMassBalance(\n", + " \tgdir_hef,\n", + " \tmb_model_class=massbalance.MonthlyTIModel,\n", + " \tmelt_f=float(melt_f),\n", + "\t\tprcp_fac=2.5,\n", + "\t\ttemp_bias=0,\n", + " \tcheck_calib_params=False,\n", + " \t) \n", + "\t\n", + "\t# Run this with the our 2 fixed mass balance parameters and our varying melt factor\n", + "\ttasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", + "\t\t\t\t\t\trun_task=tasks.run_from_climate_data, # running from observed climate data\n", + "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", + "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", + "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + "\t\t\t\t\t\tmb_model=mb_model, # use the mass balance model with the new melt_f\n", + "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", + "\t\n", + "\t# Reading the glaciological and hydrological outputs\n", + "\twith xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + "\t\t# The last step of hydrological output is NaN (we can't compute it for this year)\n", + "\t\tds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['mf_tb_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", - "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_mf_tb:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_mf_tb = ds_mf_tb.isel(time=slice(0, -1)).load()\n", - "\n", - "# Plot the runoff again for the calibrated melt_f parameter\n", - "sel_vars = [v for v in ds_mf_tb.variables if 'month_2d' not in ds_mf_tb[v].dims]\n", - "df_annual_mf_tb = ds_mf_tb[sel_vars].to_dataframe()\n", - "\n", - "mf_tb_calib" + "\t# Calculating the runoff\n", + "\tmelt_f_dict[melt_f] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier']) * 1e-9" ] }, { "cell_type": "markdown", - "id": "246378e3", + "id": "916cea1c", "metadata": {}, "source": [ - "Now calibrating `melt_f`, `prcp_fac` and `temp_bias`, we firstly calibrate `melt_f` and fix the remaining two parameters to the default values. Since we know from previous calibrations, with the values chosen for our calibration, a solution can be found with just calibrating `melt_f` alone and therefore `prcp_fac` and `temp_bias` do not reach calibration. Therefore returning the same results as before in this tutorial where `melt_f` is calibrated first:" - ] - }, - { - "cell_type": "markdown", - "id": "d58ddcbf", - "metadata": {}, - "source": [ - "### Calibrating all three parameters:\n", - "\n", - "Now we will use `mb_calibration_from_scalar_mb` to calibrate the three parameters. This works as the following:\n", - "1. Set calibrate_param1, calibrate_param2 and calibrate_param3.\n", - "2. Fix the parameters assigned to calibrate_param2 and calibrate_param3.\n", - "3. Calibrate calibrate_param1.\n", - "4. If a solution is found, stop calibration here.\n", - "5. If no solution can be found, move onto calibrate the next parameter.\n", - "6. Repeat steps 4 and 5, therefore either reaching a calibrated solution or, if no solution can be found for any parameter, an error is raised and calibration stops with no solution.\n", - "\n", - "Here we can see that the chosen order to calibrate parameters can affect the outcomes of the calibration." + "Now we plot:" ] }, { "cell_type": "code", "execution_count": null, - "id": "7077bbef", + "id": "67a19e04", "metadata": {}, "outputs": [], "source": [ - "# Calibrate all of the parameters\n", - "calib_all = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " calibrate_param1='melt_f',\n", - " calibrate_param2='prcp_fac',\n", - " calibrate_param3='temp_bias',\n", - " overwrite_gdir=True)\n", - "\n", - "file_id = '_calib_all'\n", + "# let's get a nice colormap centered around melt_f\n", + "norm = matplotlib.colors.Normalize(vmin=1.5, vmax=17)\n", + "colors_melt_f = plt.get_cmap('coolwarm')\n", "\n", - "# Adding the mass balance model with the new calibration\n", - "mbmod = massbalance.MonthlyTIModel(gdir_hef)\n", - "fls = gdir_hef.read_pickle('inversion_flowlines')\n", - "mbdf['calib_all_mb'] = mbmod.get_specific_mb(fls=fls, year=mbdf.index)\n", + "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", + "for j, melt_f in enumerate(melt_f_dict.keys()):\n", + " axs[0].plot(melt_f, melt_f_dict[melt_f].mean(), 'o',\n", + " color=colors_melt_f(norm(melt_f))) \n", + "axs[0].set_ylabel('Runoff (Mt/yr)')\n", + "axs[0].set_xlabel('Melt factor');\n", + "axs[0].set_title('Mean Runoff vs Melt Factor')\n", "\n", - "# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "# Run this again with the calibrated parameters\n", - "tasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - " run_task=tasks.run_from_climate_data, # running from observed climate data\n", - " ys=2000, # Period which we will average and constantly repeat\n", - " init_model_yr=2000, # Start from spinup year 2000\n", - " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", - " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", "\n", - "with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_calib:\n", - " # The last step of hydrological output is NaN (we can't compute it for this year)\n", - " ds_calib = ds_calib.isel(time=slice(0, -1)).load()\n", + "for melt_f in melt_f_dict.keys():\n", + " axs[1].plot(np.arange(2000,2020,1),\n", + " melt_f_dict[melt_f].values, '-', \n", + " color=colors_melt_f(norm(melt_f)),\n", + " label=melt_f)\n", + "axs[1].set_ylabel('Runoff (Mt/yr)')\n", + "axs[1].set_xlabel('Year')\n", + "axs[1].legend(title='Melt factor:', bbox_to_anchor=(1,1))\n", + "axs[1].set_xticks(np.arange(2000,2020,2));\n", + "axs[1].set_title('Runoff Time Series for Different Melt Factors')\n", "\n", - "# Plot the runoff again for the calibrated melt_f parameter\n", - "sel_vars = [v for v in ds_calib.variables if 'month_2d' not in ds_calib[v].dims]\n", - "df_annual_calib = ds_calib[sel_vars].to_dataframe()\n", - "\n", - "calib_all" + "fig.suptitle(\"Runoff and Melt Factor - {}\".format(gdir_hef.rgi_id));\n", + "plt.tight_layout()" ] }, { "cell_type": "markdown", - "id": "b7e6bb8f", + "id": "f76c8cb8", "metadata": {}, "source": [ - "We can see above that the calibrated parameters are now the same as in cases where we calibrate the `melt_f` first, and when we calibrate `melt_f` and `prcp_fac`. This is due to the behaviour of the scalar calibration function, which calibrates one parameter at a time in the order that we pass into the calibration function.\n", - "\n", - "This also means that the order in which we calibrate the parameters will change the outcome of our calibration, so keep this in mind when using the `mb_calibration_from_scalar_mb` function. In this tutorial we have kept to the same order of calibration in the permutation:\n", - "1. `melt_f`\n", - "2. `prcp_fac`\n", - "3. `temp_bias`\n", - "\n", - "But this order can be changed and is up to the discretion of the user to decide what they would like to calibrate initially." + "Again, we can see that the runoff is sensitive to the melt factor too! And that, as the melt factor increases, so does the runoff." ] }, { "cell_type": "markdown", - "id": "31fa073d", + "id": "91e6f879", "metadata": {}, "source": [ - "Now lets investigate the effect of the parameter calibration on the runoff, the aim of this tutorial!\n", - "\n", - "First we will convert the runoff values to megatonnes for readability." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0288677b", - "metadata": {}, - "outputs": [], - "source": [ - "# Convert them to megatonnes (instead of kg)\n", - "df_runoff_melt_f = df_annual_melt_f[runoff_vars] * 1e-9\n", - "df_runoff_prcp_fac = df_annual_prcp_fac[runoff_vars] * 1e-9\n", - "df_runoff_temp_bias = df_annual_temp_bias[runoff_vars] * 1e-9\n", - "df_calib = df_annual_calib[runoff_vars] * 1e-9\n", - "df_runoff_mf_pf = df_annual_mf_pf[runoff_vars] * 1e-9\n", - "df_runoff_pf_tb = df_annual_pf_tb[runoff_vars] * 1e-9\n", - "df_runoff_mf_tb = df_annual_mf_tb[runoff_vars] * 1e-9" + "We have explored the relationship between the mass balance parameters and the runoff, but what about other glaciohydrological outputs?" ] }, { "cell_type": "markdown", - "id": "1c34f50d", - "metadata": {}, - "source": [ - "And plotting runoff output values resulting from all calibration of parameters..." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "87e51a4a", + "id": "fdbca36b", "metadata": {}, - "outputs": [], "source": [ - "fig, ax = plt.subplots(figsize=(10, 3.5), sharex=True)\n", - "df_calib.sum(axis=1).plot(ax=ax, label='all params calibrated');\n", - "df_runoff_melt_f.sum(axis=1).plot(ax=ax, label='calibrated melt_f');\n", - "df_runoff.sum(axis=1).plot(ax=ax, label='default calib');\n", - "df_runoff_prcp_fac.sum(axis=1).plot(ax=ax, label='calibrated prcp_fac');\n", - "df_runoff_temp_bias.sum(axis=1).plot(ax=ax, label='calibrated temp_bias');\n", - "df_runoff_mf_pf.sum(axis=1).plot(ax=ax, label='melt_f & prcp_fac calibrated', linestyle='dashed');\n", - "df_runoff_pf_tb.sum(axis=1).plot(ax=ax, label='prcp_fac & temp_bias calibrated', linestyle='dashed');\n", - "df_runoff_mf_tb.sum(axis=1).plot(ax=ax, label='melt_f & temp_bias calibrated', linestyle='dashed');\n", - "plt.ylabel('Runoff in Megatonnes'); plt.xlabel('Years'); plt.title(f'Total annual runoff for {rgi_id}');\n", - "plt.legend();" + "# Exploring other Glaciohydrological outputs in the OGGM" ] }, { "cell_type": "markdown", - "id": "5955278f", + "id": "398c38af", "metadata": {}, "source": [ - "The impact of calibrating different permutations of parameters can be seen above, this seems to impact the runoff result quite significantly depending on the calibration choices made! It can be seen that the rough shape is similar for all calibrations, but that different permutation of calibrating parameters can have a large impact on the amount of runoff predicted by the OGGM.\n", + "**First let's start with a definition!** In this next section we will be investigating the melt contribution to runoff, defined as: $ \\frac{\\text{melt on glacier}}{\\text{runoff}}$. \n", "\n", - "We can see that some of these runoff values are the same for certain calibrations, note that here that this is due to the style of calibration of our parameters, as we can see above the calibrations that calibrate `prcp_fac` first have a significantly lower runoff than the rest." + "This represents how much of the total runoff can be attributed to the glacial melt (i.e. melt water produced from the currently glaciated area)" ] }, { "cell_type": "markdown", - "id": "8bad322a", + "id": "fe5a38fe", "metadata": {}, "source": [ - "Now plotting the mass-balance results together and comparing this to the runoff outputs:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c9455a05", - "metadata": {}, - "outputs": [], - "source": [ - "# Now plotting all the mass balance results together\n", + "Below we will investigate both the melt contribution, and `melt_on_glacier` sensitivity to the precipitation factor. \n", "\n", - "# Adding the in-situ mass balance observations for comparison\n", - "mbdf['in_situ_mb'] = gdir_hef.get_ref_mb_data().loc[2000:2019]['ANNUAL_BALANCE']\n", - "\n", - "# The mass-balance from diffrent calibration runs throughout the tutorial\n", - "plt.plot(mbdf['in_situ_mb'], label='In-situ MB')\n", - "plt.plot(mbdf['melt_f_mb'], label='melt_f')\n", - "plt.plot(mbdf['prcp_fac_mb'], label='prcp_fac')\n", - "plt.plot(mbdf['temp_bias_mb'], label='temp_bias')\n", - "plt.plot(mbdf['calib_all_mb'], label='All params', linestyle='dashed')\n", - "plt.plot(mbdf['mf_pf_mb'], label='melt_f & prcp_fac', linestyle='dashed')\n", - "plt.plot(mbdf['pf_tb_mb'], label='prcp_fac & temp_bias', linestyle='dashed')\n", - "plt.plot(mbdf['mf_tb_mb'], label='melt_f & temp_bias', linestyle='dashed')\n", - "plt.legend()\n", - "plt.xlabel(\"Year\")\n", - "plt.ylabel(\"Mass Balance\");\n", - "plt.title(f'Mass Balance time series for {rgi_id}');" - ] - }, - { - "cell_type": "markdown", - "id": "e4aa15ef", - "metadata": {}, - "source": [ - "In this graph we can now see that our different calibrations also have an impact on the mass-balance\n", - "\n", - "Here, we can see that there are a few overlaps again for different permutations of calibrations, if we look closely we can see that these are the same as the above for runoff.\n", - "\n", - "However, from initial inspection we can see that for the runoff there is quite a significant difference based off of choices of calibrations. Here we can see that the differences between the mass-balance don't seem quite as dramatic.\n", - "\n", - "From the runoff graph, the calibration of the `prcp_fac` appears to be impactful on the runoff. We can now explore a bit further, and investigate the relationship between `prcp_fac` and `melt_on_glacier`, which is a component of the `runoff`. We will do this by fixing the `prcp_fac` to a range of values and calibrating the `melt_f` to obtain a range of calibrated values. This will allow us to verify the sensitivity of these parameters." - ] - }, - { - "cell_type": "markdown", - "id": "eef65f7e", - "metadata": {}, - "source": [ - "## Sensitivity Analysis of runoff parameters in calibration\n", - "\n", - "**Sensitivity Analysis** is a technique used to determine how a model output is affected by the changes in parameters. We will perform a simple sensitvity analysis in this tutorial to investigate the affects of parameter calibration on some of the runoff components.\n", - "\n", - "We will now investigate the sensitivity of the runoff parameters to the calibration of the mass-balance parameters below. This will allow us to understand the relationship between parameters further, for example if one parameter is changed, how much does this affect the other parameter?" - ] - }, - { - "cell_type": "markdown", - "id": "af1a9529", - "metadata": {}, - "source": [ - "**First let's start with some definitions!**\n", - "\n", - "We will be investigating the melt contribution, which is:\n", - "$ \\frac{\\text{melt on glacier}}{\\text{runoff`}}$" - ] - }, - { - "cell_type": "markdown", - "id": "003110aa", - "metadata": {}, - "source": [ - "In this tutorial we will investigate the sensitivities of the `melt_on_glacier` and melt contribution to the `prcp_fac`.\n", - "\n", - "We have chosen these parameters as `melt_on_glacier` is a single variable, defined by a linear relationship from the `prcp_fac`. Here we will just be affecting the `melt_on_glacier` directly.\n", - "\n", - "The melt contribution variable is a bit more complicated, as this will be investigating the relationship between the run off which includes `melt_on_glacier` variable." + "We run a final set of simulations to obtain the glaciohydrological outputs from the OGGM for our sensitivity study." ] }, { "cell_type": "code", "execution_count": null, - "id": "834456bf", + "id": "ddd3e653", "metadata": {}, "outputs": [], "source": [ "# Varying prcp_fac between a range of values with a step of 1.0\n", - "pd_prcp_sens = pd.DataFrame(index=np.arange(0.1, 5.0, 0.3))\n", + "pd_prcp_sens = pd.DataFrame(index=np.arange(0.1, 10.0, 0.5))\n", "file_id = '_sens'\n", "\n", - "# Calibrate the melt factor for each precipitation factor\n", - "spec_mb_melt_f_sens_dict = {}\n", "for pf in pd_prcp_sens.index:\n", - " calib_param = mb_calibration_from_scalar_mb(gdir_hef,\n", - " ref_mb = ref_mb,\n", - " prcp_fac = pf,\n", - " ref_period=cfg.PARAMS['geodetic_mb_period'],\n", - " overwrite_gdir=True)\n", "\n", - " # Fill the dataframe with the calibrated parameters\n", - " pd_prcp_sens.loc[pf, 'melt_f'] = calib_param['melt_f']\n", + " mb_model = MultipleFlowlineMassBalance(\n", + " gdir_hef,\n", + " mb_model_class=massbalance.MonthlyTIModel,\n", + " prcp_fac=float(pf),\n", + " melt_f=5.0,\n", + " temp_bias=0,\n", + " check_calib_params=False,\n", + " ) \n", "\n", " # We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", " # Run this again with the calibrated parameters\n", @@ -826,6 +508,7 @@ " ys=2000, # Period which we will average and constantly repeat\n", " init_model_yr=2000, # Start from spinup year 2000\n", " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", + " mb_model=mb_model, # use the mass balance model with the new melt_f\n", " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", "\n", @@ -837,123 +520,140 @@ " sel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", " df_annual_sens = ds_sens[sel_vars].to_dataframe()\n", "\n", - " pd_prcp_sens.loc[pf, 'melt_off_glacier'] = df_annual_sens['melt_off_glacier'].mean()\n", - " pd_prcp_sens.loc[pf, 'melt_on_glacier'] = df_annual_sens['melt_on_glacier'].mean()\n", - " pd_prcp_sens.loc[pf, 'liq_prcp_off_glacier'] = df_annual_sens['liq_prcp_off_glacier'].mean()\n", - " pd_prcp_sens.loc[pf, 'liq_prcp_on_glacier'] = df_annual_sens['liq_prcp_on_glacier'].mean()\n", - " pd_prcp_sens.loc[pf, 'runoff'] = pd_prcp_sens.loc[pf, 'melt_off_glacier'] + pd_prcp_sens.loc[pf, 'melt_on_glacier'] + pd_prcp_sens.loc[pf, 'liq_prcp_off_glacier'] + pd_prcp_sens.loc[pf, 'liq_prcp_on_glacier']\n", - "\n", - "colors_melt_f = plt.get_cmap('viridis').colors[10::10]\n", - "plt.figure()\n", - "for j, pf in enumerate(pd_prcp_sens.index):\n", - " plt.plot(pf, pd_prcp_sens.loc[pf, 'melt_on_glacier'], 'o', color=colors_melt_f[j])\n", - " plt.xlabel('Precipitation factor')\n", - " plt.ylabel('Mean annual melt on glacier (kg m-2 yr-1)')\n", - " plt.title('Sensitivity of Melt On Glacier to Precipitation Factor')\n", - "\n", - "# Here we are varying the precipitation factor and calibrating the melt factor for each value\n", - "# Then plotting the mean annual melt_on_glacier against the precipitation factor" + " pd_prcp_sens.loc[pf, 'melt_off_glacier'] = df_annual_sens['melt_off_glacier'].mean() * 1e-9\n", + " pd_prcp_sens.loc[pf, 'melt_on_glacier'] = df_annual_sens['melt_on_glacier'].mean() * 1e-9\n", + " pd_prcp_sens.loc[pf, 'liq_prcp_off_glacier'] = df_annual_sens['liq_prcp_off_glacier'].mean() * 1e-9\n", + " pd_prcp_sens.loc[pf, 'liq_prcp_on_glacier'] = df_annual_sens['liq_prcp_on_glacier'].mean() * 1e-9\n", + " pd_prcp_sens.loc[pf, 'runoff'] = pd_prcp_sens.loc[pf, 'melt_off_glacier'] + pd_prcp_sens.loc[pf, 'melt_on_glacier'] + pd_prcp_sens.loc[pf, 'liq_prcp_off_glacier'] + pd_prcp_sens.loc[pf, 'liq_prcp_on_glacier']" ] }, { "cell_type": "markdown", - "id": "00bed96c", + "id": "611d9aa2", "metadata": {}, "source": [ - "In the above plot, a strong linear positive correlation is shown between the precipitation factor and the mean annual melt on the glacier. \n", + "Now plotting to investigate the sensitivity of the `melt_on_glacier` and the melt contribution to the precipitation factor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8650f395", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", + "\n", + "for pf in pd_prcp_sens.index:\n", + " axs[0].scatter(\n", + " pf,\n", + " pd_prcp_sens.loc[pf, 'melt_on_glacier'], # melt on glacier\n", + " color='blue'\n", + " )\n", + " \n", + " axs[1].scatter(\n", + " pf,\n", + " pd_prcp_sens.loc[pf, 'melt_on_glacier']/pd_prcp_sens.loc[pf, 'runoff'], # melt contribution\n", + " color='blue'\n", + " )\n", "\n", - "From the calibration, it can be seen that as the precipitation factor is changed, the melt factor parameter is calibrated accordingly and increases with the precipitation factor to maintain the average annual mass-balance. Therefore as the `prcp_fac` increases, the `melt_on_glacier` increases linearly in a 1:1 ratio.\n", "\n", - "In the OGGM, the `melt_on_glacier` is derived from both the `melt_f` parameter and `prcp_fac`. Therefore this sensitivity makes sense." + "axs[0].set_xlabel('Precipitation factor')\n", + "axs[0].set_ylabel('Mean annual melt on glacier (Mt/yr)')\n", + "axs[0].set_title('Sensitivity of Melt on Glacier to Precipitation Factor')\n", + "\n", + "axs[1].set_xlabel('Precipitation factor')\n", + "axs[1].set_ylabel('Mean Glacial Melt Contribution')\n", + "axs[1].set_title('Sensitivity of Glacial Melt Contribution to the Precipitation Factor')\n", + "\n", + "plt.show()\n" ] }, { "cell_type": "markdown", - "id": "28875fe5", + "id": "ceb52f26", "metadata": {}, "source": [ - "Now investigating the sensitivity of the Glacier Melt Contribution to the Precipitation Factor." + "We can see that as the precipitation factor increases, the mean annual melt on glacier increases. \n", + "\n", + "However, when we divide this value by the total runoff to derive the melt contribution, this results in a decrease in the melt contribution as the precipitation factor increases.\n", + "\n", + "Let's think about why there might be an increase in the melt on glacier with an increasing precipitation factor, but a decreasing glacial melt contribution? We can start by investigating how much the total runoff increases, compared to the melt on glacier, and this might help us answer the question!\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "7ad8fa07", + "id": "e33f3221", "metadata": {}, "outputs": [], "source": [ - "colors_melt_f = plt.get_cmap('viridis').colors[10::10]\n", - "plt.figure()\n", - "for j, pf in enumerate(pd_prcp_sens.index):\n", - " plt.plot(pf, pd_prcp_sens.loc[pf, 'melt_on_glacier']/pd_prcp_sens.loc[pf, 'runoff'], 'o', color=colors_melt_f[j])\n", - " plt.xlabel('Precipitation factor')\n", - " plt.ylabel('Mean Glacial Melt Contribution')\n", - " plt.title('Sensitivity of Glacial Melt Contribution to the Precipitation Factor')\n", + "fig, axs = plt.subplots()\n", + "\n", + "for i, pf in enumerate(pd_prcp_sens.index):\n", + " axs.scatter(\n", + " pf,\n", + " pd_prcp_sens.loc[pf, 'melt_on_glacier'],\n", + " color='blue',\n", + " label='Melt on glacier' if i == 0 else None\n", + " )\n", + " axs.scatter(\n", + " pf,\n", + " pd_prcp_sens.loc[pf, 'runoff'],\n", + " color='red',\n", + " label='Runoff' if i == 0 else None\n", + " )\n", + " \n", + "axs.set_xlabel('Precipitation factor')\n", + "axs.set_ylabel('Glaciohydrological components')\n", + "axs.set_title('Sensitivity of Melt on Glacier to Precipitation Factor')\n", + "\n", + "plt.legend()\n", "\n", - "# Here we are plotting how the glacial melt contribution (melt_on_glacier/runoff) varies with the changing precipitation factor\n", - "# This shows how the relative contribution of glacier melt to total runoff changes as we adjust precipitation" + "plt.show()\n" ] }, { "cell_type": "markdown", - "id": "75444665", + "id": "f1594fc4", "metadata": {}, "source": [ - "It can be seen that there is a relatively strong negative correlation, between the precipitation factor and the mean glacial melt contribution. And therefore it can be seen that the mean glacial melt contribution is sensitive to the change in precipitation factor. As we increase the `prcp_fac`, we know that the `melt_f` increases along with this. Considering the melt contribution variable contains a numerator of `melt_on_glacier`, which has been shown to increase as the `prcp_fac` parameter increases.\n", - "\n", - "The denominator is the runoff, which contains 4 summed variables (including `melt_on_glacier`) that are affected by the `prcp_fac`, all of which increase as the `prcp_fac` increases, as can be seen within the OGGM model implementation. The runoff therefore increases at least as much as the `melt_on_glacier`.\n", + "We can see that the total runoff is more sensitive to the changes in the precipitation factor and is increasing at a much more rapid rate when we increase the precipitation factor. Therefore this decrease in the melt contribution makes sense.\n", "\n", - "Therefore the increase of the denominator overpowers the increase of the numerator, and the relationship shown above reflects this reaction to the increase in `prcp_fac`." + "This shows us that different outputs can have different sensitivities to the same changing input parameters, and can lead to some interesting results!" ] }, { "cell_type": "markdown", - "id": "cf62afed", + "id": "dcdb8013", "metadata": {}, "source": [ - "## Exercise:\n", + "# Tutorial take-aways\n", + "- Sensitivity analysis is a useful tool to understand the relationship between our model inputs and outputs\n", + "- One-at-a-time sensitivity analysis is an exploratory tool to help us better investigate model behaviour.\n", + "- Different outputs can exhibit different sensitivities for the same range of input values.\n", + "- The model sensitivity depends greatly on the chosen inputs (what they are, and their ranges) and the outputs we are considering.\n", "\n", - "Try this yourself for sensitivity with the precipitation factor when varying the melt factor as above:\n", - "1. The melt_off_glacier parameter.\n", - "2. The total liquid precipitation parameters (sum of on and off).\n", - "3. Any other interesting combinations you can discover!\n", "\n", - "Does anything you discover surprise you?\n", - "Are there particularly sensitive parameters compared to others?\n", - "\n", - "**Stretch:** Can you experiment different combinations of parameters for sensitvity, e.g. such as varying the temp_bias parameters instead? " + "This is a very simple application of sensitivity analysis applied to the OGGM, but these insights can motivate further use of sensitivity analysis on a larger-scale. There are future plans to integrate the [Sensitivity Analysis For Everyone Toolbox (SAFE)](https://safetoolbox.github.io/) with OGGM for more complex sensitivity analysis applications!\n" ] }, { "cell_type": "markdown", - "id": "8e0eae95", + "id": "c2c538a5", "metadata": {}, "source": [ - "## WGMS" + "# References" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "53b9a503", + "cell_type": "markdown", + "id": "37cd2b95", "metadata": {}, - "outputs": [], "source": [ - "# # We start by fetching the observations\n", - "# mbdf = gdir_hef.get_ref_mb_data()[['ANNUAL_BALANCE']]\n", - "# mbdf.columns = ['in_situ_mb']\n", - "# mbdf = mbdf.loc[2000:2019]\n", - "# mbdf.plot();\n", - "# mbdf.mean()" + "1. Francesca Pianosi, Keith Beven, Jim Freer, Jim W. Hall, Jonathan Rougier, David B. Stephenson, Thorsten Wagener,\n", + "Sensitivity analysis of environmental models: A systematic review with practical workflow, Environmental Modelling & Software, Volume 79, 2016, Pages 214-232, ISSN 1364-8152, https://doi.org/10.1016/j.envsoft.2016.02.008." ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "155dbd5a", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { From 15a6e81d4014d7556afa8b38f7cbcb756682cb60 Mon Sep 17 00:00:00 2001 From: Chloe Hancock Date: Fri, 10 Apr 2026 09:32:26 +0100 Subject: [PATCH 4/5] Fixed typos/spelling --- notebooks/tutorials/runoff_sensitivity.ipynb | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/notebooks/tutorials/runoff_sensitivity.ipynb b/notebooks/tutorials/runoff_sensitivity.ipynb index edb788e5..98fe23f5 100644 --- a/notebooks/tutorials/runoff_sensitivity.ipynb +++ b/notebooks/tutorials/runoff_sensitivity.ipynb @@ -158,7 +158,7 @@ "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly outputs provide additional information\n", "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", "\t\n", "\t# Now read the hydrological outputs for this run\n", @@ -169,7 +169,7 @@ "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", "\n", - "\t# Calcuating the runoff for the temperature bias value and storing it in a dictionary to plot\n", + "\t# Calculating the runoff for the temperature bias value and storing it in a dictionary to plot\n", "\ttemp_bias_dict[temp_bias] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier'])*1e-9 # Convert from Kilograms to Megatonnes per year (Mt/yr) for easier plotting" ] }, @@ -237,7 +237,7 @@ "id": "8ad1e197", "metadata": {}, "source": [ - "Now lets explore what happens when we very the **precipitation factor**, and fix our other two mass balance parameters." + "Now lets explore what happens when we vary the **precipitation factor**, and fix our other two mass balance parameters." ] }, { @@ -269,7 +269,7 @@ "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly outputs provide additional information\n", "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", "\t\n", "\t# Reading the glaciological and hydrological outputs for this run again\n", @@ -377,7 +377,7 @@ "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", "\t\t\t\t\t\tmb_model=mb_model, # use the mass balance model with the new melt_f\n", - "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly outputs provide additional information\n", "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", "\t\n", "\t# Reading the glaciological and hydrological outputs\n", @@ -466,7 +466,7 @@ "source": [ "**First let's start with a definition!** In this next section we will be investigating the melt contribution to runoff, defined as: $ \\frac{\\text{melt on glacier}}{\\text{runoff}}$. \n", "\n", - "This represents how much of the total runoff can be attributed to the glacial melt (i.e. melt water produced from the currently glaciated area)" + "This represents how much of the total runoff can be attributed to the glacial melt (i.e. the melt water produced from the currently glaciated area)." ] }, { @@ -509,7 +509,7 @@ " init_model_yr=2000, # Start from spinup year 2000\n", " init_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", " mb_model=mb_model, # use the mass balance model with the new melt_f\n", - " store_monthly_hydro=True, # Monthly ouptuts provide additional information\n", + " store_monthly_hydro=True, # Monthly outputs provide additional information\n", " output_filesuffix=file_id); # an identifier for the output file, to read it later\n", "\n", " with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", @@ -574,6 +574,8 @@ "id": "ceb52f26", "metadata": {}, "source": [ + "The above graph shows that these glaciohydrological outputs are both sensitive to the precipitation factor changing!\n", + "\n", "We can see that as the precipitation factor increases, the mean annual melt on glacier increases. \n", "\n", "However, when we divide this value by the total runoff to derive the melt contribution, this results in a decrease in the melt contribution as the precipitation factor increases.\n", @@ -618,7 +620,7 @@ "id": "f1594fc4", "metadata": {}, "source": [ - "We can see that the total runoff is more sensitive to the changes in the precipitation factor and is increasing at a much more rapid rate when we increase the precipitation factor. Therefore this decrease in the melt contribution makes sense.\n", + "We can see that the total runoff is more sensitive to the changes in the precipitation factor and is increasing at a much more rapid rate when we increase the precipitation factor. Therefore this causes a decrease in the melt contribution makes sense.\n", "\n", "This shows us that different outputs can have different sensitivities to the same changing input parameters, and can lead to some interesting results!" ] From a4c79b64c958fb1ef41fdb722928c30164917e17 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Fri, 10 Apr 2026 23:03:41 +0100 Subject: [PATCH 5/5] cosmetic changed --- notebooks/tutorials/runoff_sensitivity.ipynb | 478 +++++++++++-------- 1 file changed, 268 insertions(+), 210 deletions(-) diff --git a/notebooks/tutorials/runoff_sensitivity.ipynb b/notebooks/tutorials/runoff_sensitivity.ipynb index 98fe23f5..1d85d12d 100644 --- a/notebooks/tutorials/runoff_sensitivity.ipynb +++ b/notebooks/tutorials/runoff_sensitivity.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "349a2484", + "id": "0", "metadata": {}, "source": [ "# Runoff Sensitivity to Mass Balance Parameters" @@ -10,38 +10,42 @@ }, { "cell_type": "markdown", - "id": "0f9dec41", + "id": "1", "metadata": {}, "source": [ - "Goals of this notebook: To explore the sensitivity of the runoff and OGGM's hydrological components to the mass balance parameters." + "Goals of this notebook: To explore the sensitivity of the runoff and OGGM's hydrological components to the mass balance model parameters.\n", + "\n", + "Mass balance parameters directly control the glacier's mass balance (accumulation and ablation), but their impact extends beyond solid ice and snow. In OGGM, runoff is computed from multiple components: glacier melt, precipitation both on and off the glacier, and snow melt. Because runoff is a sum of these components (and some components are more sensitive to mass balance parameter changes than others), **runoff can exhibit even greater sensitivity to parameter variations than the mass balance itself**.\n", + "\n", + "Recent work ([Wimberly et al., 2025](https://tc.copernicus.org/articles/19/1491/2025/)) has shown that these sensitivities have important implications for glacier hydrology and projections. Understanding how mass balance parameters propagate through the hydrological cycle is critical for water resource assessments and future projections.\n", + "\n", + "This will allow us to understand the relationship between our parameters and model output, and to appreciate why careful parameter calibration is important for accurate hydrological predictions." ] }, { "cell_type": "markdown", - "id": "543e58a1", + "id": "2", "metadata": {}, "source": [ - "Our other hydrology tutorials introduce the idea of hydrology in the OGGM. It is recommended to consult the previous tutorials to gain an understanding first as to how hydrology works in the OGGM. \n", - "\n", - "[Glaciers as water resources: Part 1](https://edu-notebooks.oggm.org/oggm-edu/glacier_water_resources.html)\n", - "\n", - "[Glaciers as water resources: Part 2](https://edu-notebooks.oggm.org/oggm-edu/glacier_water_resources_projections.html)" + "It is recommended to consult the previous tutorials to gain an understanding first as to how runoff is computed in OGGM:\n", + "- [Glaciers as water resources: Part 1](https://edu-notebooks.oggm.org/oggm-edu/glacier_water_resources.html)\n", + "- [Glaciers as water resources: Part 2](https://edu-notebooks.oggm.org/oggm-edu/glacier_water_resources_projections.html)" ] }, { "cell_type": "markdown", - "id": "d1abcbfd", + "id": "3", "metadata": {}, "source": [ "## Set Up\n", "\n", - "First install required packages to run this tutorial and initialize our glacier directories!" + "First import the required packages to run this tutorial and initialize our glacier directories!" ] }, { "cell_type": "code", "execution_count": null, - "id": "0b6e8193", + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -59,7 +63,7 @@ { "cell_type": "code", "execution_count": null, - "id": "85804d49", + "id": "5", "metadata": {}, "outputs": [], "source": [ @@ -70,43 +74,45 @@ }, { "cell_type": "markdown", - "id": "9d2efb09", + "id": "6", "metadata": {}, "source": [ - "We start from a well known glacier in the Austrian Alps, Hintereisferner. But you can choose any other glacier, e.g. from [this list](https://github.com/OGGM/oggm-sample-data/blob/master/wgms/rgi_wgms_links_20220112.csv)" + "We start from a well known glacier in the Austrian Alps, Hintereisferner. But you can choose any other glacier, e.g. from [this list](https://github.com/OGGM/oggm-sample-data/blob/master/wgms/rgi_wgms_links_20220112.csv)." ] }, { "cell_type": "code", "execution_count": null, - "id": "402e44d9", + "id": "7", "metadata": {}, "outputs": [], "source": [ "# Hintereisferner\n", "rgi_id = 'RGI60-11.00897'\n", "\n", - "# We pick the elevation-bands glaciers because they run a bit faster - but they create more step changes in the area outputs\n", - "gdir_hef = workflow.init_glacier_directories([rgi_id], from_prepro_level=4, prepro_border=160, prepro_base_url=DEFAULT_BASE_URL)[0]" + "# We pick the elevation-bands glaciers because they run a bit faster -\n", + "# but they create more step changes in the area outputs\n", + "gdir_hef = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=160,\n", + " prepro_base_url=DEFAULT_BASE_URL)[0]" ] }, { "cell_type": "markdown", - "id": "eef65f7e", + "id": "8", "metadata": {}, "source": [ "# An Introduction to Sensitivity Analysis\n", "\n", - "**Sensitivity Analysis** investigates how the variation in the output of a numerical model can be attributed to variations of its input factors [1](https://www.sciencedirect.com/science/article/pii/S1364815216300287).\n", + "**Sensitivity Analysis** investigates how the variation in the output of a numerical model can be attributed to variations of its input factors ([Pianosi et al., (2016)](https://www.sciencedirect.com/science/article/pii/S1364815216300287)).\n", "\n", "In this tutorial, we perform a simple, exploratory one-at-a-time sensitivity analysis to investigate the effects of the mass balance parameters on the glaciohydrological model outputs. We will focus on the mass balance parameters: the melt factor, temperature bias and precipitation factor.\n", "\n", - "This will allow us to understand the relationship between our parameters and model output, for example if one parameter is changed, how much does this affect our model output?" + "**Important: unlike an \"operational\" OGGM run, the parameters below are varied without ensuring that observations are matched. This would lead to different results / conclusions and could be the topic of another tutorial.**" ] }, { "cell_type": "markdown", - "id": "bb462087", + "id": "9", "metadata": {}, "source": [ "# Simple Sensitivity Analysis Experiment" @@ -114,7 +120,7 @@ }, { "cell_type": "markdown", - "id": "050f0271", + "id": "10", "metadata": {}, "source": [ "### Temperature Bias" @@ -122,7 +128,7 @@ }, { "cell_type": "markdown", - "id": "a9c0ce89", + "id": "11", "metadata": {}, "source": [ "We will begin by investigating the sensitvity of the runoff to one parameter at a time, we will start with **temperature bias**. Below, we will vary only the temperature bias, and fix the melt factor and the precipitation factor." @@ -131,51 +137,61 @@ { "cell_type": "code", "execution_count": null, - "id": "2d1cf29c", + "id": "12", "metadata": {}, "outputs": [], "source": [ - "temp_bias_dict = {}\n", + "temp_bias_df = pd.DataFrame()\n", "file_id = '_sens'\n", "\n", - "for temp_bias in np.arange(-5,5.0,0.5): # We are varying the temperature bias\n", - "\n", - "\t# For each temperature bias, we create a mass balance model with the same melt factor and precipitation factor, but with the new temperature bias\n", - "\tmb_model = MultipleFlowlineMassBalance(\n", - " \tgdir_hef,\n", - " \tmb_model_class=massbalance.MonthlyTIModel,\n", - " \ttemp_bias=float(temp_bias), # Vary the temperature bias\n", - "\t\tmelt_f=5.0, # Fix melt factor to 5.0 for all runs, to see the effect of the temperature bias alone\n", - "\t\tprcp_fac=2.5, # Fix precipitation factor to 2.5 for all runs, to see the effect of the temperature bias alone\n", - " \tcheck_calib_params=False,\n", - " \t) \n", - "\t\n", - "\t# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "\t# Run this with the our 2 fixed mass balance parameters and our varying temperature bias\n", - "\ttasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - "\t\t\t\t \t\trun_task=tasks.run_from_climate_data, # Run from climate data\n", - "\t\t\t\t \t\tmb_model=mb_model, # use the mass balance model with the new temp_bias\n", - "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", - "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", - "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly outputs provide additional information\n", - "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", - "\t\n", - "\t# Now read the hydrological outputs for this run\n", - "\twith xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", - "\t\t# The last step of hydrological output is NaN (we can't compute it for this year)\n", - "\t\tds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", - "\n", - "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", - "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", - "\n", - "\t# Calculating the runoff for the temperature bias value and storing it in a dictionary to plot\n", - "\ttemp_bias_dict[temp_bias] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier'])*1e-9 # Convert from Kilograms to Megatonnes per year (Mt/yr) for easier plotting" + "for temp_bias in np.arange(-5, 5.0, 0.5): # We are varying the temperature bias\n", + "\n", + " # For each temperature bias, we create a mass balance model with the same melt factor\n", + " # and precipitation factor, but with the new temperature bias\n", + " mb_model = MultipleFlowlineMassBalance(\n", + " gdir_hef,\n", + " mb_model_class=massbalance.MonthlyTIModel,\n", + " temp_bias=float(temp_bias), # Vary the temperature bias\n", + " melt_f=5.0, # Fix melt factor to 5.0 for all runs\n", + " prcp_fac=2.5, # Fix precipitation factor to 2.5 for all runs\n", + " check_calib_params=False, # We are forcing parameters which do not match observations\n", + " )\n", + "\n", + " # We are using the task run_with_hydro to store hydrological outputs along with\n", + " # the usual glaciological outputs\n", + " tasks.run_with_hydro(\n", + " gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # Run from climate data\n", + " mb_model=mb_model, # Use the mass balance model with the new temp_bias\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # Use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly outputs provide additional information\n", + " output_filesuffix=file_id, # Identifier for the output file, to read it later\n", + " )\n", + "\n", + " # Now read the hydrological outputs for this run\n", + " with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", + "\n", + " sel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + " df_annual_sens = ds_sens[sel_vars].to_dataframe()\n", + "\n", + " # Store the runoff time series in a DataFrame, one column per temperature bias\n", + " temp_bias_df[temp_bias] = (\n", + " df_annual_sens['melt_off_glacier']\n", + " + df_annual_sens['melt_on_glacier']\n", + " + df_annual_sens['liq_prcp_off_glacier']\n", + " + df_annual_sens['liq_prcp_on_glacier']\n", + " ) * 1e-9 # Convert from kilograms to megatonnes per year (Mt/yr) for easier plotting\n", + "\n", + "temp_bias_df = temp_bias_df.sort_index(axis=1)" ] }, { "cell_type": "markdown", - "id": "96c586d6", + "id": "13", "metadata": {}, "source": [ "Now, for each temperature bias value, let's plot the mean runoff against the temperature bias, and plot the runoff time series to undestand how this varies across our range of temperature bias values." @@ -184,7 +200,7 @@ { "cell_type": "code", "execution_count": null, - "id": "801bc112", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -192,41 +208,47 @@ "norm = matplotlib.colors.Normalize(vmin=-5, vmax=5.01)\n", "colors_temp_bias = plt.get_cmap('coolwarm')\n", "\n", - "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", - "for j, temp_bias in enumerate(temp_bias_dict.keys()):\n", - " axs[0].plot(temp_bias, temp_bias_dict[temp_bias].mean(), 'o',\n", - " color=colors_temp_bias(norm(temp_bias))) \n", + "fig, axs = plt.subplots(1, 2, figsize=(14, 6))\n", + "for temp_bias in temp_bias_df.columns:\n", + " axs[0].plot(\n", + " temp_bias,\n", + " temp_bias_df[temp_bias].mean(),\n", + " 'o',\n", + " color=colors_temp_bias(norm(temp_bias)),\n", + " )\n", "axs[0].set_ylabel('Mean Runoff (Mt/yr)')\n", - "axs[0].set_xlabel('temp_bias (°C)');\n", + "axs[0].set_xlabel('temp_bias (°C)')\n", "axs[0].set_title('Mean Runoff vs Temperature Bias')\n", "\n", - "\n", - "for temp_bias in temp_bias_dict.keys():\n", - " axs[1].plot(np.arange(2000,2020,1),\n", - " temp_bias_dict[temp_bias].values, '-', \n", - " color=colors_temp_bias(norm(temp_bias)),\n", - " label=temp_bias)\n", + "for temp_bias in temp_bias_df.columns:\n", + " axs[1].plot(\n", + " temp_bias_df.index,\n", + " temp_bias_df[temp_bias].values,\n", + " '-',\n", + " color=colors_temp_bias(norm(temp_bias)),\n", + " label=temp_bias,\n", + " )\n", "axs[1].set_ylabel('Runoff (Mt/yr)')\n", "axs[1].set_xlabel('Year')\n", - "axs[1].legend(title='temp_bias:', bbox_to_anchor=(1,1))\n", - "axs[1].set_xticks(np.arange(2000,2020,2));\n", + "axs[1].legend(title='temp_bias:', bbox_to_anchor=(1, 1))\n", + "axs[1].set_xticks(temp_bias_df.index[::2])\n", "axs[1].set_title('Runoff Time Series for Different Temperature Biases')\n", "\n", - "fig.suptitle(\"Runoff and Temperature Bias - {}\".format(gdir_hef.rgi_id));\n", + "fig.suptitle(f'Runoff and Temperature Bias - {gdir_hef.rgi_id}')\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", - "id": "e6917d8d", + "id": "15", "metadata": {}, "source": [ - "We can see that the runoff appears to be sensitive to the temperature bias! In the graph on the left, we can see the mean runoff increases as we increase the temperature bias, and on the graph on the right, we can see how this affects the runoff annually." + "We can see that the runoff appears to be sensitive to the temperature bias! In the graph on the left, we can see the mean runoff increases as we increase the temperature bias (likely because of an increase of glacier melt), and on the graph on the right, we can see how this affects the runoff annually. The high temperature bias leads to a strong reduction of runoff towards the end of the period, likely because the glacier is melting very fast and cannot sustain high runoff for long." ] }, { "cell_type": "markdown", - "id": "e2461a95", + "id": "16", "metadata": {}, "source": [ "### Precipitation Factor" @@ -234,7 +256,7 @@ }, { "cell_type": "markdown", - "id": "8ad1e197", + "id": "17", "metadata": {}, "source": [ "Now lets explore what happens when we vary the **precipitation factor**, and fix our other two mass balance parameters." @@ -243,49 +265,60 @@ { "cell_type": "code", "execution_count": null, - "id": "b659f889", + "id": "18", "metadata": {}, "outputs": [], "source": [ - "prcp_fac_dict = {}\n", - "\n", - "for prcp_fac in np.arange(0.1,10,0.5): # Now we are varying the precipitation factor\n", - "\t\n", - "\t# For each precipitation factor, we create a mass balance model with the same melt factor and temperature bias, but with the new precipitation factor\n", - "\tmb_model = MultipleFlowlineMassBalance(\n", - " \tgdir_hef,\n", - " \tmb_model_class=massbalance.MonthlyTIModel,\n", - " \tprcp_fac=float(prcp_fac),\n", - "\t\tmelt_f=5.0, # Fix melt factor\n", - "\t\ttemp_bias=0, # Fix the temperature bias to 0\n", - " \tcheck_calib_params=False,\n", - " \t) \n", - "\t\n", - "\t# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", - "\t# Run this with the our 2 fixed mass balance parameters and our varying precipitation factor\n", - "\ttasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - "\t\t\t\t \t\trun_task=tasks.run_from_climate_data,\n", - "\t\t\t\t \t\tmb_model=mb_model, # use the mass balance model with the new prcp_fac\n", - "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", - "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", - "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly outputs provide additional information\n", - "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", - "\t\n", - "\t# Reading the glaciological and hydrological outputs for this run again\n", - "\twith xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", - "\t\t# The last step of hydrological output is NaN (we can't compute it for this year)\n", - "\t\tds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", - "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", - "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", - "\n", - "\t# Calculating the runoff\n", - "\tprcp_fac_dict[prcp_fac] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier']) * 1e-9" + "prcp_fac_df = pd.DataFrame()\n", + "\n", + "for prcp_fac in np.arange(0.1, 10, 0.5): # Now we are varying the precipitation factor\n", + "\n", + " # For each precipitation factor, we create a mass balance model with the same melt\n", + " # factor and temperature bias, but with the new precipitation factor\n", + " mb_model = MultipleFlowlineMassBalance(\n", + " gdir_hef,\n", + " mb_model_class=massbalance.MonthlyTIModel,\n", + " prcp_fac=float(prcp_fac),\n", + " melt_f=5.0, # Fix melt factor\n", + " temp_bias=0, # Fix the temperature bias to 0\n", + " check_calib_params=False, # We are forcing parameters which do not match observations\n", + " )\n", + "\n", + " # We are using the task run_with_hydro to store hydrological outputs along with\n", + " # the usual glaciological outputs\n", + " tasks.run_with_hydro(\n", + " gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data,\n", + " mb_model=mb_model, # Use the mass balance model with the new prcp_fac\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # Use the previous run as initial state\n", + " store_monthly_hydro=True, # Monthly outputs provide additional information\n", + " output_filesuffix=file_id, # Identifier for the output file, to read it later\n", + " )\n", + "\n", + " # Reading the glaciological and hydrological outputs for this run again\n", + " with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", + "\n", + " sel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + " df_annual_sens = ds_sens[sel_vars].to_dataframe()\n", + "\n", + " # Store the runoff time series in a DataFrame, one column per precipitation factor\n", + " prcp_fac_df[prcp_fac] = (\n", + " df_annual_sens['melt_off_glacier']\n", + " + df_annual_sens['melt_on_glacier']\n", + " + df_annual_sens['liq_prcp_off_glacier']\n", + " + df_annual_sens['liq_prcp_on_glacier']\n", + " ) * 1e-9\n", + "\n", + "prcp_fac_df = prcp_fac_df.sort_index(axis=1)" ] }, { "cell_type": "markdown", - "id": "ae390ad3", + "id": "19", "metadata": {}, "source": [ "Now let's plot to see our results." @@ -294,7 +327,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9ff119cd", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -302,40 +335,47 @@ "norm = matplotlib.colors.Normalize(vmin=0.1, vmax=10)\n", "colors_prcp_fac = plt.get_cmap('coolwarm')\n", "\n", - "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", - "for j, prcp_fac in enumerate(prcp_fac_dict.keys()):\n", - " axs[0].plot(prcp_fac, prcp_fac_dict[prcp_fac].mean(), 'o',\n", - " color=colors_prcp_fac(norm(prcp_fac))) \n", + "fig, axs = plt.subplots(1, 2, figsize=(14, 6))\n", + "for prcp_fac in prcp_fac_df.columns:\n", + " axs[0].plot(\n", + " prcp_fac,\n", + " prcp_fac_df[prcp_fac].mean(),\n", + " 'o',\n", + " color=colors_prcp_fac(norm(prcp_fac)),\n", + " )\n", "axs[0].set_ylabel('Mean Runoff (Mt/yr)')\n", - "axs[0].set_xlabel('Precipitation factor');\n", + "axs[0].set_xlabel('Precipitation factor')\n", "axs[0].set_title('Mean Runoff vs Precipitation Factor')\n", "\n", - "for prcp_fac in prcp_fac_dict.keys():\n", - " axs[1].plot(np.arange(2000,2020,1),\n", - " prcp_fac_dict[prcp_fac].values, '-', \n", - " color=colors_prcp_fac(norm(prcp_fac)),\n", - " label=prcp_fac)\n", + "for prcp_fac in prcp_fac_df.columns:\n", + " axs[1].plot(\n", + " prcp_fac_df.index,\n", + " prcp_fac_df[prcp_fac].values,\n", + " '-',\n", + " color=colors_prcp_fac(norm(prcp_fac)),\n", + " label=prcp_fac,\n", + " )\n", "axs[1].set_ylabel('Runoff (Mt/yr)')\n", "axs[1].set_xlabel('Year')\n", - "axs[1].legend(title='Precipitation factor:', bbox_to_anchor=(1,1))\n", - "axs[1].set_xticks(np.arange(2000,2020,2));\n", + "axs[1].legend(title='Precipitation factor:', bbox_to_anchor=(1, 1))\n", + "axs[1].set_xticks(prcp_fac_df.index[::2])\n", "axs[1].set_title('Runoff Time Series for Different Precipitation Factors')\n", "\n", - "fig.suptitle(\"Runoff and Precipitation Factor - {}\".format(gdir_hef.rgi_id));\n", + "fig.suptitle(f'Runoff and Precipitation Factor - {gdir_hef.rgi_id}')\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", - "id": "44c2ecc9", + "id": "21", "metadata": {}, "source": [ - "We can see that again, the runoff appears sensitive to the precipitation factor! And that as the precipitation factor increases, as does the runoff." + "We can see that again, the runoff appears sensitive to the precipitation factor! And that as the precipitation factor increases, as does the runoff. Note that the increase is not linear." ] }, { "cell_type": "markdown", - "id": "0958f452", + "id": "22", "metadata": {}, "source": [ "### Melt Factor" @@ -343,7 +383,7 @@ }, { "cell_type": "markdown", - "id": "2764c928", + "id": "23", "metadata": {}, "source": [ "Finally, let's see what happens when we alter the melt factor and fix the remaining 2 mass balance parameters!" @@ -352,49 +392,59 @@ { "cell_type": "code", "execution_count": null, - "id": "eb55224d", + "id": "24", "metadata": {}, "outputs": [], "source": [ - "melt_f_dict = {}\n", - "\n", - "for melt_f in np.arange(1.5,17,1.0): # Now we are varying the melt factor\n", - "\t\n", - "\t# For each melt factor, we create a mass balance model with the same precipitation factor and temperature bias, but now with the new melt factor\n", - "\tmb_model = MultipleFlowlineMassBalance(\n", - " \tgdir_hef,\n", - " \tmb_model_class=massbalance.MonthlyTIModel,\n", - " \tmelt_f=float(melt_f),\n", - "\t\tprcp_fac=2.5,\n", - "\t\ttemp_bias=0,\n", - " \tcheck_calib_params=False,\n", - " \t) \n", - "\t\n", - "\t# Run this with the our 2 fixed mass balance parameters and our varying melt factor\n", - "\ttasks.run_with_hydro(gdir_hef, # Run on the selected glacier\n", - "\t\t\t\t\t\trun_task=tasks.run_from_climate_data, # running from observed climate data\n", - "\t\t\t\t\t\tys=2000, # Period which we will average and constantly repeat\n", - "\t\t\t\t\t\tinit_model_yr=2000, # Start from spinup year 2000\n", - "\t\t\t\t\t\tinit_model_filesuffix='_spinup_historical', # use the previous run as initial state\n", - "\t\t\t\t\t\tmb_model=mb_model, # use the mass balance model with the new melt_f\n", - "\t\t\t\t\t\tstore_monthly_hydro=True, # Monthly outputs provide additional information\n", - "\t\t\t\t\t\toutput_filesuffix=file_id) # an identifier for the output file, to read it later\n", - "\t\n", - "\t# Reading the glaciological and hydrological outputs\n", - "\twith xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", - "\t\t# The last step of hydrological output is NaN (we can't compute it for this year)\n", - "\t\tds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", - "\n", - "\tsel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", - "\tdf_annual_sens = ds_sens[sel_vars].to_dataframe()\n", - "\n", - "\t# Calculating the runoff\n", - "\tmelt_f_dict[melt_f] = (df_annual_sens['melt_off_glacier']+df_annual_sens['melt_on_glacier']+df_annual_sens['liq_prcp_off_glacier']+ df_annual_sens['liq_prcp_on_glacier']) * 1e-9" + "melt_f_df = pd.DataFrame()\n", + "\n", + "for melt_f in np.arange(1.5, 17, 1.0): # Now we are varying the melt factor\n", + "\n", + " # For each melt factor, we create a mass balance model with the same precipitation\n", + " # factor and temperature bias, but now with the new melt factor\n", + " mb_model = MultipleFlowlineMassBalance(\n", + " gdir_hef,\n", + " mb_model_class=massbalance.MonthlyTIModel,\n", + " melt_f=float(melt_f),\n", + " prcp_fac=2.5,\n", + " temp_bias=0,\n", + " check_calib_params=False, # We are forcing parameters which do not match observations\n", + " )\n", + "\n", + " # Run this with our 2 fixed mass balance parameters and our varying melt factor\n", + " tasks.run_with_hydro(\n", + " gdir_hef, # Run on the selected glacier\n", + " run_task=tasks.run_from_climate_data, # Running from observed climate data\n", + " ys=2000, # Period which we will average and constantly repeat\n", + " init_model_yr=2000, # Start from spinup year 2000\n", + " init_model_filesuffix='_spinup_historical', # Use the previous run as initial state\n", + " mb_model=mb_model, # Use the mass balance model with the new melt_f\n", + " store_monthly_hydro=True, # Monthly outputs provide additional information\n", + " output_filesuffix=file_id, # Identifier for the output file, to read it later\n", + " )\n", + "\n", + " # Reading the glaciological and hydrological outputs\n", + " with xr.open_dataset(gdir_hef.get_filepath('model_diagnostics', filesuffix=file_id)) as ds_sens:\n", + " # The last step of hydrological output is NaN (we can't compute it for this year)\n", + " ds_sens = ds_sens.isel(time=slice(0, -1)).load()\n", + "\n", + " sel_vars = [v for v in ds_sens.variables if 'month_2d' not in ds_sens[v].dims]\n", + " df_annual_sens = ds_sens[sel_vars].to_dataframe()\n", + "\n", + " # Store the runoff time series in a DataFrame, one column per melt factor\n", + " melt_f_df[melt_f] = (\n", + " df_annual_sens['melt_off_glacier']\n", + " + df_annual_sens['melt_on_glacier']\n", + " + df_annual_sens['liq_prcp_off_glacier']\n", + " + df_annual_sens['liq_prcp_on_glacier']\n", + " ) * 1e-9\n", + "\n", + "melt_f_df = melt_f_df.sort_index(axis=1)" ] }, { "cell_type": "markdown", - "id": "916cea1c", + "id": "25", "metadata": {}, "source": [ "Now we plot:" @@ -403,7 +453,7 @@ { "cell_type": "code", "execution_count": null, - "id": "67a19e04", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -411,41 +461,47 @@ "norm = matplotlib.colors.Normalize(vmin=1.5, vmax=17)\n", "colors_melt_f = plt.get_cmap('coolwarm')\n", "\n", - "fig, axs = plt.subplots(1,2,figsize=(14,6))\n", - "for j, melt_f in enumerate(melt_f_dict.keys()):\n", - " axs[0].plot(melt_f, melt_f_dict[melt_f].mean(), 'o',\n", - " color=colors_melt_f(norm(melt_f))) \n", + "fig, axs = plt.subplots(1, 2, figsize=(14, 6))\n", + "for melt_f in melt_f_df.columns:\n", + " axs[0].plot(\n", + " melt_f,\n", + " melt_f_df[melt_f].mean(),\n", + " 'o',\n", + " color=colors_melt_f(norm(melt_f)),\n", + " )\n", "axs[0].set_ylabel('Runoff (Mt/yr)')\n", - "axs[0].set_xlabel('Melt factor');\n", + "axs[0].set_xlabel('Melt factor')\n", "axs[0].set_title('Mean Runoff vs Melt Factor')\n", "\n", - "\n", - "for melt_f in melt_f_dict.keys():\n", - " axs[1].plot(np.arange(2000,2020,1),\n", - " melt_f_dict[melt_f].values, '-', \n", - " color=colors_melt_f(norm(melt_f)),\n", - " label=melt_f)\n", + "for melt_f in melt_f_df.columns:\n", + " axs[1].plot(\n", + " melt_f_df.index,\n", + " melt_f_df[melt_f].values,\n", + " '-',\n", + " color=colors_melt_f(norm(melt_f)),\n", + " label=melt_f,\n", + " )\n", "axs[1].set_ylabel('Runoff (Mt/yr)')\n", "axs[1].set_xlabel('Year')\n", - "axs[1].legend(title='Melt factor:', bbox_to_anchor=(1,1))\n", - "axs[1].set_xticks(np.arange(2000,2020,2));\n", + "axs[1].legend(title='Melt factor:', bbox_to_anchor=(1, 1))\n", + "axs[1].set_xticks(melt_f_df.index[::2])\n", "axs[1].set_title('Runoff Time Series for Different Melt Factors')\n", "\n", - "fig.suptitle(\"Runoff and Melt Factor - {}\".format(gdir_hef.rgi_id));\n", + "fig.suptitle(f'Runoff and Melt Factor - {gdir_hef.rgi_id}')\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", - "id": "f76c8cb8", + "id": "27", "metadata": {}, "source": [ - "Again, we can see that the runoff is sensitive to the melt factor too! And that, as the melt factor increases, so does the runoff." + "Again, we can see that the runoff is sensitive to the melt factor too! And that, as the melt factor increases, so does the runoff. Like with the temperature bias, high melt factors lead to a quick depletion of the glacier." ] }, { "cell_type": "markdown", - "id": "91e6f879", + "id": "28", "metadata": {}, "source": [ "We have explored the relationship between the mass balance parameters and the runoff, but what about other glaciohydrological outputs?" @@ -453,40 +509,42 @@ }, { "cell_type": "markdown", - "id": "fdbca36b", + "id": "29", "metadata": {}, "source": [ - "# Exploring other Glaciohydrological outputs in the OGGM" + "# Exploring other Glaciohydrological outputs in OGGM" ] }, { "cell_type": "markdown", - "id": "398c38af", + "id": "30", "metadata": {}, "source": [ - "**First let's start with a definition!** In this next section we will be investigating the melt contribution to runoff, defined as: $ \\frac{\\text{melt on glacier}}{\\text{runoff}}$. \n", + "**First let's start with a definition!** In this next section we will be investigating the melt contribution to runoff, defined as: \n", + "\n", + "$ \\frac{\\text{melt on glacier}}{\\text{runoff}}$. \n", "\n", "This represents how much of the total runoff can be attributed to the glacial melt (i.e. the melt water produced from the currently glaciated area)." ] }, { "cell_type": "markdown", - "id": "fe5a38fe", + "id": "31", "metadata": {}, "source": [ "Below we will investigate both the melt contribution, and `melt_on_glacier` sensitivity to the precipitation factor. \n", "\n", - "We run a final set of simulations to obtain the glaciohydrological outputs from the OGGM for our sensitivity study." + "We run a final set of simulations to obtain the glaciohydrological outputs from OGGM for our sensitivity study." ] }, { "cell_type": "code", "execution_count": null, - "id": "ddd3e653", + "id": "32", "metadata": {}, "outputs": [], "source": [ - "# Varying prcp_fac between a range of values with a step of 1.0\n", + "# Varying prcp_fac between a range of values with a step of 0.5\n", "pd_prcp_sens = pd.DataFrame(index=np.arange(0.1, 10.0, 0.5))\n", "file_id = '_sens'\n", "\n", @@ -499,7 +557,7 @@ " melt_f=5.0,\n", " temp_bias=0,\n", " check_calib_params=False,\n", - " ) \n", + " )\n", "\n", " # We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs\n", " # Run this again with the calibrated parameters\n", @@ -529,7 +587,7 @@ }, { "cell_type": "markdown", - "id": "611d9aa2", + "id": "33", "metadata": {}, "source": [ "Now plotting to investigate the sensitivity of the `melt_on_glacier` and the melt contribution to the precipitation factor." @@ -538,7 +596,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8650f395", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -550,7 +608,7 @@ " pd_prcp_sens.loc[pf, 'melt_on_glacier'], # melt on glacier\n", " color='blue'\n", " )\n", - " \n", + "\n", " axs[1].scatter(\n", " pf,\n", " pd_prcp_sens.loc[pf, 'melt_on_glacier']/pd_prcp_sens.loc[pf, 'runoff'], # melt contribution\n", @@ -571,7 +629,7 @@ }, { "cell_type": "markdown", - "id": "ceb52f26", + "id": "35", "metadata": {}, "source": [ "The above graph shows that these glaciohydrological outputs are both sensitive to the precipitation factor changing!\n", @@ -586,7 +644,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e33f3221", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -605,7 +663,7 @@ " color='red',\n", " label='Runoff' if i == 0 else None\n", " )\n", - " \n", + "\n", "axs.set_xlabel('Precipitation factor')\n", "axs.set_ylabel('Glaciohydrological components')\n", "axs.set_title('Sensitivity of Melt on Glacier to Precipitation Factor')\n", @@ -617,7 +675,7 @@ }, { "cell_type": "markdown", - "id": "f1594fc4", + "id": "37", "metadata": {}, "source": [ "We can see that the total runoff is more sensitive to the changes in the precipitation factor and is increasing at a much more rapid rate when we increase the precipitation factor. Therefore this causes a decrease in the melt contribution makes sense.\n", @@ -627,7 +685,7 @@ }, { "cell_type": "markdown", - "id": "dcdb8013", + "id": "38", "metadata": {}, "source": [ "# Tutorial take-aways\n", @@ -637,12 +695,12 @@ "- The model sensitivity depends greatly on the chosen inputs (what they are, and their ranges) and the outputs we are considering.\n", "\n", "\n", - "This is a very simple application of sensitivity analysis applied to the OGGM, but these insights can motivate further use of sensitivity analysis on a larger-scale. There are future plans to integrate the [Sensitivity Analysis For Everyone Toolbox (SAFE)](https://safetoolbox.github.io/) with OGGM for more complex sensitivity analysis applications!\n" + "This is a very simple application of sensitivity analysis applied to OGGM, but these insights can motivate further use of sensitivity analysis on a larger-scale. There are future plans to integrate the [Sensitivity Analysis For Everyone Toolbox (SAFE)](https://safetoolbox.github.io/) with OGGM for more complex sensitivity analysis applications!\n" ] }, { "cell_type": "markdown", - "id": "c2c538a5", + "id": "39", "metadata": {}, "source": [ "# References" @@ -650,11 +708,11 @@ }, { "cell_type": "markdown", - "id": "37cd2b95", + "id": "40", "metadata": {}, "source": [ - "1. Francesca Pianosi, Keith Beven, Jim Freer, Jim W. Hall, Jonathan Rougier, David B. Stephenson, Thorsten Wagener,\n", - "Sensitivity analysis of environmental models: A systematic review with practical workflow, Environmental Modelling & Software, Volume 79, 2016, Pages 214-232, ISSN 1364-8152, https://doi.org/10.1016/j.envsoft.2016.02.008." + "- Pianosi, F., Beven, K., Freer, J., Hall, J. W., Rougier, J., Stephenson, D. B., and Wagener, T.: Sensitivity analysis of environmental models: a systematic review with practical workflow, Environ. Model. Softw., 79, 214–232, https://doi.org/10.1016/j.envsoft.2016.02.008, 2016\n", + "- Wimberly, F., Ultee, L., Schuster, L., Huss, M., Rounce, D. R., Maussion, F., Coats, S., Mackay, J., and Holmgren, E.: Inter-model differences in 21st century glacier runoff for the world's major river basins, The Cryosphere, 19, 1491–1511, https://doi.org/10.5194/tc-19-1491-2025, 2025. " ] } ], @@ -669,7 +727,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.14" + "version": "3.14.0" } }, "nbformat": 4,