diff --git a/_toc.yml b/_toc.yml index 174698a..712a20e 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 0000000..0184422 --- /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 4a2399d..61a349a 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)" ] }, {