diff --git a/notebooks/tutorials/massbalance_calibration.ipynb b/notebooks/tutorials/massbalance_calibration.ipynb index 2376910..aa0501d 100644 --- a/notebooks/tutorials/massbalance_calibration.ipynb +++ b/notebooks/tutorials/massbalance_calibration.ipynb @@ -13,7 +13,7 @@ "source": [ "The default mass-balance (MB) model of OGGM is a very standard [temperature index melt model](https://www.sciencedirect.com/science/article/pii/S0022169403002579). \n", "\n", - "In versions before 1.6, OGGM had a complex calibration procedure which originated from the times where we had only observations from a few hundred glaciers. We used them to calibrate the model and then a so-called infamous *tstar* (infamous in very niche circles of first- and second-gen OGGMers) which was interpolated to glaciers without observations (see the [original publication](https://www.the-cryosphere.net/6/1295/2012/tc-6-1295-2012.html)). This method was very powerful but, as new observational datasets emerged, we can now calibrate on a glacier-per-glacier basis. With the new era of geodetic observations, OGGM uses the average geodetic observations from Jan 2000--Jan 2020 of [Hugonnet al. 2021](https://www.nature.com/articles/s41586-021-03436-z), that are now available for almost every glacier world-wide. \n", + "In versions before 1.6, OGGM had a complex calibration procedure which originated from the times, where we had only observations from a few hundred glaciers. We used them to calibrate the model and then a so-called infamous *tstar* (infamous in very niche circles of first- and second-gen OGGMers) which was interpolated to glaciers without observations (see the [original publication](https://www.the-cryosphere.net/6/1295/2012/tc-6-1295-2012.html)). This method was very powerful but, as new observational datasets emerged, we can now calibrate on a glacier-per-glacier basis. With the new era of geodetic observations, OGGM uses the average geodetic observations from Jan 2000--Jan 2020 of [Hugonnet al. 2021](https://www.nature.com/articles/s41586-021-03436-z), that are now available for almost every glacier world-wide.\n", "\n", "Pre-processed directories from OGGM (from the Bremen server) have been calibrated for you, based on a specific climate dataset (W5E5) and our own dedicated calibration strategy. But, what if you want to use another climate dataset? Or another reference dataset?\n", "\n", @@ -40,9 +40,7 @@ "import matplotlib\n", "import pandas as pd\n", "import numpy as np\n", - "import os\n", "\n", - "import oggm\n", "from oggm import cfg, utils, workflow, tasks, graphics\n", "from oggm.core import massbalance\n", "from oggm.core.massbalance import mb_calibration_from_scalar_mb, mb_calibration_from_geodetic_mb, mb_calibration_from_wgms_mb" @@ -56,14 +54,14 @@ "source": [ "cfg.initialize(logging_level='WARNING')\n", "cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-calib-mb', reset=True)\n", - "cfg.PARAMS['border'] = 10" + "cfg.PARAMS['border'] = 80" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We start from two well known glaciers in the Austrian Alps, Kesselwandferner and Hintereisferner. But you can also choose any other 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 two well known glaciers in the Austrian Alps, Kesselwandferner and Hintereisferner. But you can also choose any other glacier, e.g. from [this list](https://github.com/OGGM/oggm-sample-data/blob/master/wgms/rgi_wgms_links_20220112.csv)." ] }, { @@ -73,7 +71,7 @@ "outputs": [], "source": [ "# we start from preprocessing level 3\n", - "base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5/'\n", + "base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2025.6/elev_bands/W5E5/per_glacier'\n", "gdirs = workflow.init_glacier_directories(['RGI60-11.00787', 'RGI60-11.00897'], from_prepro_level=3, prepro_base_url=base_url)" ] }, @@ -137,7 +135,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This is called the \"baseline climate\" in OGGM and is necessary to calibrate the model against observations. Ideally, the baseline climate should be real observations as perfect as possible, but in reality this is not the case. Often, gridded climate datasets have biases - we need to take this into accound during our calibration. Let's have a look at the mass balance parameters for both glaciers:" + "This is called the \"baseline climate\" in OGGM and is necessary to calibrate the model against observations. Ideally, the baseline climate should be real observations as perfect as possible, but in reality this is not the case. Often, gridded climate datasets have biases - we need to take this into account during our calibration. Let's have a look at the mass balance parameters for both glaciers:" ] }, { @@ -166,11 +164,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will explain later what all these mean. Lets focus on these for now: `melt_f`, `prcp_fac`, and `temp_bias` which have been calibrated with `reference_mb` over the reference period `reference_period`.\n", + "We will explain later what all these mean. Let's focus on these for now: `melt_f`, `prcp_fac`, and `temp_bias` which have been calibrated with `reference_mb` over the reference period `reference_period`.\n", "\n", "Per default the [Hugonnet et al. (2021)](https://www.nature.com/articles/s41586-021-03436-z) average geodetic observation is used over the entire time period Jan 2000 to Jan 2020 to calibrate the MB model parameter(s) for every single glacier.\n", "\n", - "For our two neighboring example glaciers, that share the same climate gridpoint, the same pre-calibrated `temp_bias` and the same `melt_f` is used. The `prcp_fac` is slighly different, hence in that example, changing the `prcp_fac` within the narrow range of $[0.8,1.2]$ was sufficient to match the MB model to the observations. \n", + "For our two neighboring example glaciers, that share the same climate gridpoint, the same pre-calibrated `temp_bias` and the same `melt_f` is used. The `prcp_fac` is slightly different, hence in that example, changing the `prcp_fac` within the narrow range of $[0.8,1.2]$ was sufficient to match the MB model to the observations.\n", "\n", "Note that the two glaciers are in the same climate (from the forcing data) but are very different in size, orientation and geometry. So, at least one of the MB model parameters is different (here the `prcp_fac`) and also the observed MB values vary:" ] @@ -212,7 +210,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note alst that there are some global MB parameters (`mb_global_params`), which we assume to be the same globally for every glacier:" + "Note also that there are some global MB parameters (`mb_global_params`), which we assume to be the same globally for every glacier:" ] }, { @@ -295,7 +293,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Perfect! Our MB model reproduces the average observed MB, so the calibrated worked!" + "Perfect! Our MB model reproduces the average observed MB, so the calibration worked!" ] }, { @@ -323,7 +321,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For example, let's calibrate the `melt_f` to match the same observations as before (i.e., average geodetic observation) and fix the other parameters. We will use the more general function [mb_calibration_from_scalar_mb](https://docs.oggm.org/en/latest/generated/oggm.tasks.mb_calibration_from_scalar_mb.html) which is most flexible. However, it requires you to give manually the observed MB on which you want to calibrate:" + "For example, let's calibrate the `melt_f` to match the same observations as before (i.e., average geodetic observation) and fix the other parameters. We will use the more general function [mb_calibration_from_scalar_mb](https://docs.oggm.org/en/latest/generated/oggm.tasks.mb_calibration_from_scalar_mb.html) which is more flexible. However, it requires you to give manually the observed MB on which you want to calibrate:" ] }, { @@ -379,7 +377,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Only `melt_f` is used to calibrate the MB model. `prcp_fac` was chosen depending on the average winter precipitation (explained below) and `temp_bias` is 0. " + "Only `melt_f` is used to calibrate the MB model. `prcp_fac` was chosen depending on the average winter precipitation (explained below) and `temp_bias` is 0." ] }, { @@ -537,7 +535,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Ok, but actually you might trust the in-situ observations much more than the geodetic observation. So if you are only interested in glaciers with these in-situ observations, you can also use the in-situ observations for calibration. This might allow you also to calibrate over a longer time period and to use additional informations (such as the interannual MB variability, seasonal MB, ...). However, bare in mind that there are often gaps in the in-situ MB time series and make sure that you calibrate to the same modelled time period!\n", + "Ok, but actually you might trust the in-situ observations much more than the geodetic observation. So if you are only interested in glaciers with these in-situ observations, you can also use the in-situ observations for calibration. This might allow you also to calibrate over a longer time period and to use additional information (such as the interannual MB variability, seasonal MB, ...). However, bear in mind that there are often gaps in the in-situ MB time series and make sure that you calibrate to the same modelled time period!\n", "\n", "Attention: For the Hintereisferner glacier with an almost 70-year long time series, the assumption that we make, i.e., that the area does not change over the time period, gets more problematic. So think twice before repeating this at home!" ] @@ -715,7 +713,7 @@ "source": [ "If we believe that our observations are true, maybe the climate or the processes represented by the MB model are erroneous. \n", "\n", - "What can we do to still match the `ref_mb`? We can either change the `melt_f` ranges by setting other values to the parameters above or we can allow that another MB model parameter is changed (i.e., `calibrate_param2`). This is basically very similar to the three-step-calibration \n", + "What can we do to still match the `ref_mb`? We can either change the `melt_f` ranges by setting other values to the parameters above, or we can allow that another MB model parameter is changed (i.e., `calibrate_param2`). This is basically very similar to the three-step-calibration\n", "first introduced in [Huss & Hock 2015](https://doi.org/10.3389/feart.2015.00054), but here you can choose your parameter ranges and parameter order yourself. \n", "\n", "To reduce these MB model calibration errors, you can first change the `melt_f`, fix it at the lower or upper limit, and then change the `temp_bias`." @@ -753,7 +751,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the same way, you can also use your own MB observations or fake data to calibrate the MB model. We use here as an example, an unrealistically high positive MB for the Hintereisferner over a 10-year time period. To allow the calibration to happen, we need again to set `calibrate_param2`! Here we will use the `prcp_fac` as second parameter for a change:" + "In the same way, you can also use your own MB observations or fake data to calibrate the MB model. We use here as an example, an unrealistically high positive MB for the Hintereisferner over a 10-year period. To allow the calibration to happen, we need again to set `calibrate_param2`! Here we will use the `prcp_fac` as second parameter for a change:" ] }, { @@ -763,7 +761,7 @@ "outputs": [], "source": [ "ref_period = '2000-01-01_2010-01-01'\n", - "ref_mb = 2000 # Let's use an unrealistically positive mass-balance\n", + "ref_mb = 2000 # Let's use an unrealistically positive mass-balance\n", "mb_calibration_from_scalar_mb(gdir_hef, ref_mb=ref_mb,\n", " ref_period=ref_period, \n", " overwrite_gdir=True, \n", @@ -871,14 +869,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Overparameteristion or the magic choice of the best calibration option:" + "## Overparameterisation or the magic choice of the best calibration option:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We found already some combinations that equally well match the average MB over a given time period. As we only use only one observation per glacier (i.e., per default the average geodetic MB from 2000-2020), but have up to three free MB model parameters, the MB model is overparameterised. That means, there are in theory an infinite amount of calibration options possible that equally well match the one obervation. Let's look a bit more systematically into that:\n", + "We found already some combinations that equally well match the average MB over a given time period. As we only use one observation per glacier (i.e., per default the average geodetic MB from 2000-2020), but have up to three free MB model parameters, the MB model is overparameterised. That means, there are in theory an infinite amount of calibration options possible that equally well match the one observation. Let's look a bit more systematically into that:\n", "\n", "We will use a range of different `prcp_fac` and then calibrate the `melt_f` accordingly to always match to the default average MB (`ref_mb`) over the reference period (`ref_period`)." ] @@ -977,7 +975,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The larger the `prcp_fac` and the `melt_f`, the larger is the interannual MB variability. For glaciers with in-situ observations, we can find a combination of `prcp_fac` and `melt_f` that has a similar interannnual MB variability than the observations. For example, you can choose a MB model parameter combination where the standard deviation quotient of the annual the modelled and observed MB is near to 1: " + "The larger the `prcp_fac` and the `melt_f`, the larger is the interannual MB variability. For glaciers with in-situ observations, we can find a combination of `prcp_fac` and `melt_f` that has a similar interannual MB variability than the observations. For example, you can choose a MB model parameter combination where the standard deviation quotient of the annual modelled and observed MB is near to 1:" ] }, { @@ -1024,7 +1022,7 @@ " mb_temp_b_sens = massbalance.MonthlyTIModel(gdir_hef)\n", " # ok, we actually matched the new ref_mb\n", " spec_mb_temp_b_sens_dict[temp_bias] = mb_temp_b_sens.get_specific_mb(h, w, year=np.arange(2000,2020,1))\n", - " except RuntimeError: #, 'RGI60-11.00897: ref mb not matched. Try to set calibrate_param2'\n", + " except RuntimeError: # 'RGI60-11.00897: ref mb not matched. Try to set calibrate_param2'\n", " pass" ] }, @@ -1096,7 +1094,7 @@ " mb_pf_temp_b_sens = massbalance.MonthlyTIModel(gdir_hef)\n", " # ok, we actually matched the new ref_mb\n", " spec_mb_pf_temp_b_sens_dict[temp_bias] = mb_pf_temp_b_sens.get_specific_mb(h, w, year=np.arange(2000,2020,1))\n", - " except RuntimeError: #, 'RGI60-11.00897: ref mb not matched. Try to set calibrate_param2'\n", + " except RuntimeError: # 'RGI60-11.00897: ref mb not matched. Try to set calibrate_param2'\n", " pass" ] }, @@ -1169,7 +1167,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the [preprocessed directories](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/) MB model calibration is done by using `mb_calibration_from_geodetic_mb` to determine the MB model parameters. We will reproduce the option from the preprocessed directory first:" + "In the [preprocessed directories](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2025.6) MB model calibration is done by using `mb_calibration_from_geodetic_mb` to determine the MB model parameters. We will reproduce the option from the preprocessed directory first:" ] }, { @@ -1222,7 +1220,7 @@ "source": [ "`prcp_fac` is chosen from a relation to the average winter precipitation. Fig.1 (below) shows the used relationship between winter precipitation and `prcp_fac`.\n", "\n", - "It was calibrated by adapting `melt_f` and `prcp_fac` to match the average geodetic and winter MB on around 100 glaciers with both informations available. The found relationship of decreasing `prcp_fac` for increasing winter precipitation makes sense, as glaciers with already a large winter precipitation should not be corrected with a large multiplicative `prcp_fac`.\n", + "It was calibrated by adapting `melt_f` and `prcp_fac` to match the average geodetic and winter MB on around 100 glaciers where both datasets are available. [Here](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/_notebooks/oggm_v16_winter_mb_match_to_prescribe_prcp_fac/figures/RGI60-13.05504_winter_mb_GSWP3_W5E5.png) is an example figure, which shows how the precipitation factor was chosen for such a glacier. The found relationship of decreasing `prcp_fac` for increasing winter precipitation makes sense, as glaciers with already a large winter precipitation should not be corrected with a large multiplicative `prcp_fac`.\n", "\n", "
\n", " \n", @@ -1246,14 +1244,21 @@ "prcp_fac_array = utils.clip_array(prcp_fac, r0, r1)\n", "plt.plot(w_prcp_array, prcp_fac_array)\n", "plt.xlabel(r'winter daily mean precipitation' +'\\n'+r'(kg m$^{-2}$ day$^{-1}$)')\n", - "plt.ylabel('precipitation factor (prcp_fac)'); plt.title('Fig. 1');" + "plt.ylabel('precipitation factor (prcp_fac)'); plt.title('Fig. 1 (only valid for W5E5 and the specific standard OGGM settings)');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Second step: data informed temperature bias " + "This relationship is only valid for W5E5 with the specific standard OGGM settings. We repeated the same approach with ERA5, and found that the precipitation factor decreases rather linearly with the winter precipitation for this dataset ([figure that shows the fits for both climate datasets and includes the underlying winter-MB-calibrated precipitation factors](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/_notebooks/oggm_v16_winter_mb_match_to_prescribe_prcp_fac/comparison_oggm_v16_GSWP3_W5E5_vs_ERA5_calib_log_fit_pf_distribution_change_monthly_cte_melt_f_minus_1.png)). The higher resolution of ERA5 compared to W5E5 may explain the different estimated relationships. Further details are provided in the [jupyter notebook used to estimate the fits](https://nbviewer.org/urls/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/_notebooks/oggm_v16_winter_mb_match_to_prescribe_prcp_fac/match_winter_mb_w5e5_era5.ipynb)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Second step: data informed temperature bias" ] }, { @@ -1262,7 +1267,7 @@ "source": [ "Using this data-informed precipitation factor, we then do a global calibration where the temperature bias (`temp_bias`) is calibrated, while the melt factor (`melt_f`) is fixed at 5 kg m-2 day-1 K-1 (default value based on [Schuster et al., 2023](https://doi.org/10.1017/aog.2023.57)). \n", "\n", - "The idea is that if many glaciers within the same grid point need a temperature bias to reach the obeserved MB, this indicates that a systematic correction is necessary (at least for this MB model in particular). In fact, we can plot the median bias required to match MB observations using this technique, which gives us the following plot:\n", + "The idea is that if many glaciers within the same grid point need a temperature bias to reach the observed MB, this indicates that a systematic correction is necessary (at least for this MB model in particular). In fact, we can plot the median bias required to match MB observations using this technique, which gives us the following plot (here for OGGM v1.6.1, but it is very similar in v1.6.3):\n", "\n", "![err](https://user-images.githubusercontent.com/10050469/224318400-ec1d8825-d7e7-4cdb-94f3-ebb95b8f7120.jpg)" ] @@ -1273,11 +1278,11 @@ "source": [ "The fact that the `temp_bias` parameter is spatially correlated (many regions are all blue or red) indicate that something in the data needs to be corrected for our model. It is this information that we use to inform the next step.\n", "\n", - "**The code we used for this step is available [here](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/calibration/1.6.1/). As explained above, we do a global run with fixed precip factor and melt factor, then store the resulting parameters in a csv file used by OGGM. The csv file can be found [here](https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/w5e5_temp_bias_v2023.4.csv).**\n", + "**The code we used for this step is available [here](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/calibration/1.6.3/). As explained above, we do a global run with a fixed precipitation factor and melt factor, then store the resulting parameters in a csv file used by OGGM. This is repeated for each climate dataset and RGI version. The csv file using W5E5 and RGI 6.2 can be found [here](https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/w5e5_rgi6_perglacier_temp_bias_v2025.6.2.csv).**\n", "\n", "
\n", " \n", - " Similar to the precip factor, we do not necessarily recommend our users to use this method in their own calibration effort. We don't think it's a bad method, but given the importance of calibration we think that diversity is important!\n", + " Similar to the precipitation factor, we do not necessarily recommend our users to use this method in their own calibration effort. We don't think it's a bad method, but given the importance of calibration we think that diversity is important!\n", " \n", "
" ] @@ -1295,9 +1300,11 @@ "source": [ "Finally, we now run [mb_calibration_from_scalar_mb](https://docs.oggm.org/en/latest/generated/oggm.tasks.mb_calibration_from_scalar_mb.html) again for each glacier, as follows:\n", "- use the first guess: `melt_f` = 5, `prcp_fac` = data-informed from step 1, `temp_bias` = data-informed from step 2\n", - "- if this doesn't match (this would be highly unlikely), allow `prcp_fac` to vary again between 0.8 and 1.2 times the original guess ($\\pm$20%). This is justified by the fact that the first guess for precipitation is also highly uncertain. If that worked, the calibration stops (33.6% of all glaciers worldwide are calibrated this way, for 41.1% of the total area).\n", - "- if the above did not work, allow `melt_f` to vary again. If that worked, the calibration stops (60.6% of all glaciers worldwide are calibrated this way, for 41.1% of the total area).\n", - "- finally, if the above did not work, allow `temp_bias` to vary again (5.9% of all glaciers worldwide are calibrated this way, for 2.2% of the total area).\n", + "- if this doesn't match (this would be highly unlikely), allow `prcp_fac` to vary again between 0.8 and 1.2 times the original guess ($\\pm$20%). This is justified by the fact that the first guess for precipitation is also highly uncertain. If that worked, the calibration stops (33.6% of all glaciers worldwide are calibrated this way, for 41.3% of the total area).\n", + "- if the above did not work, allow `melt_f` to vary again. If that worked, the calibration stops (60.4% of all glaciers worldwide are calibrated this way, for 56.4% of the total area).\n", + "- finally, if the above did not work, allow `temp_bias` to vary again (6.0% of all glaciers worldwide are calibrated this way, for 2.2% of the total area).\n", + "\n", + "*(numbers computed by the [massbalance_global_params.ipynb notebook](massbalance_global_params.ipynb) and are valid for the OGGM v1.6.3 (2025.6) W5E5 gdir)*\n", "\n", "
\n", " \n", @@ -1335,8 +1342,8 @@ "- We can use different observational data for calibration: \n", " - calibrating to geodetic observations using different time periods, to in-situ direct glaciological observations from the WGMS (if available) or to other custom MB data. \n", "- There exist different ways of calibrating to the average observational data:\n", - " - You can use the same option (`informed_threestep`) as done operationally for the [preprocessed directories (L>=3)](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/elev_bands/W5E5/) on all glaciers world-wide \n", - " - you can calibrate the `melt_f`, and have the `prcp_fac` and `temp_bias` fixed. If the calibration does not work, the `temp_bias` can be varied aswell. \n", + " - You can use the same option (`informed_threestep`) as done operationally for the [preprocessed directories (L>=3)](https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2025.6) on all glaciers world-wide\n", + " - you can calibrate the `melt_f`, and have the `prcp_fac` and `temp_bias` fixed. If the calibration does not work, the `temp_bias` can be varied as well.\n", " - you can also calibrate instead on `prcp_fac` or `temp_bias` and fix the other parameters.\n", "- However, we showed that the parameter combination choice has an influence on other estimates than the average MB (which then influences also future projections). \n", "- As user of OGGM, you might just use the calibrated MB model parameters from the preprocessed directories. Nevertheless, it is good to be aware of the overparameterisation problem. In addition, if you want to include uncertainties of the MB model calibration, you could include additional experiments that use another calibration option and create your own preprocessed directories with these options. With more available observational data or improved climate data, you might also be able to use better ways to calibrate the MB model parameters. \n" @@ -1356,11 +1363,6 @@ ], "metadata": { "hide_input": false, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -1371,7 +1373,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.11.3" }, "latex_envs": { "LaTeX_envs_menu_present": true, diff --git a/notebooks/tutorials/runoff_calibration.ipynb b/notebooks/tutorials/runoff_calibration.ipynb new file mode 100644 index 0000000..6a73b87 --- /dev/null +++ b/notebooks/tutorials/runoff_calibration.ipynb @@ -0,0 +1,1065 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Influence of mass-balance parameters on glacier runoff\n", + "\n", + "In this notebook, we will calabrate the mass-balance parameters in the OGGM, these parameters will then be used to calculate the components of the runoff and the runoff variable itself. These variables will then be investigated to understand sensitvity against the calibrated parameters.\n", + "\n", + "For more information on the calibration methods, we would recommend the tutorial on mass-balance calibration [massbalance_calibration.ipynb](https://tutorials.oggm.org/stable/notebooks/tutorials/massbalance_calibration.html). Here our focus is how these impact runoff, this work has been motivated by the Wimberly et al paper [Inter-model differences in 21st centry glacier runoff for the world's major river basins](https://tc.copernicus.org/articles/19/1491/2025/).\n", + "\n", + "In this notebook we will:\n", + "1. Use the OGGM scalar calibration method\n", + "2. Show how users can calibrate parameters and use these for runoff outputs to understand impact of parameter changes.\n", + "3. Investigate variable sensitivity to the calibrations." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "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": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import numpy as np\n", + "\n", + "import oggm\n", + "from oggm import cfg, utils, workflow, tasks\n", + "import xarray as xr\n", + "from oggm.core import massbalance\n", + "from oggm.core.massbalance import mb_calibration_from_scalar_mb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "cfg.initialize(logging_level='WARNING')\n", + "cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-calib-ro', reset=True)\n", + "cfg.PARAMS['store_model_geometry'] = True" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "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": "5", + "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", + "base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5_spinup'\n", + "gdir_hef = workflow.init_glacier_directories([rgi_id], from_prepro_level=4, prepro_border=160, prepro_base_url=base_url)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Generating the Hydrological and Glaciological outputs" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "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": "8", + "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": "9", + "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": "10", + "metadata": {}, + "source": [ + "The figure below shows the runoff for the default calibrated parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "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": "12", + "metadata": {}, + "source": [ + "This was done with the following parameter set:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "gdir_hef.read_json('mb_calib')" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## Parameter Calibration" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "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 the time of writing this tutorial October 2025) only 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": "16", + "metadata": {}, + "source": [ + "### **Calibration:** Scalar Mass Balance Calibration" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "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": "18", + "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": "19", + "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": "20", + "metadata": {}, + "source": [ + "### Single Parameter Calibration" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "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": "22", + "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": "23", + "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": "24", + "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": "25", + "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": "26", + "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": "27", + "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": "28", + "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": "29", + "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": "30", + "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": "31", + "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": "32", + "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": "33", + "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": "34", + "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": "35", + "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": "36", + "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": "37", + "metadata": {}, + "source": [ + "Now calibrating `prcp_fac` and `temp_bias`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38", + "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": "39", + "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": "40", + "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": "41", + "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": "42", + "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": "43", + "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": "44", + "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": "45", + "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": "46", + "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": "47", + "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": "48", + "metadata": {}, + "source": [ + "And plotting runoff output values resulting from all calibration of parameters..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49", + "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": "50", + "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": "51", + "metadata": {}, + "source": [ + "Now plotting the mass-balance results together and comparing this to the runoff outputs:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52", + "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": "53", + "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": "54", + "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": "55", + "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": "56", + "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": "57", + "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": "58", + "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": "59", + "metadata": {}, + "source": [ + "Now investigating the sensitivity of the Glacier Melt Contribution to the Precipitation Factor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60", + "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": "61", + "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": "62", + "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": "63", + "metadata": {}, + "source": [ + "## WGMS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64", + "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": "65", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "oggm_env", + "language": "python", + "name": "python3" + }, + "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.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}