From 0a2f64c2cdcfdbbdf14a4a2a83b0097b91f566f3 Mon Sep 17 00:00:00 2001 From: Alexander Fischer <42172336+afisc@users.noreply.github.com> Date: Thu, 22 Aug 2024 17:16:19 +0200 Subject: [PATCH 1/3] minor typo corrections (#122) Co-authored-by: afisc --- notebooks/tutorials/massbalance_calibration.ipynb | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/notebooks/tutorials/massbalance_calibration.ipynb b/notebooks/tutorials/massbalance_calibration.ipynb index 71e47530..23769108 100644 --- a/notebooks/tutorials/massbalance_calibration.ipynb +++ b/notebooks/tutorials/massbalance_calibration.ipynb @@ -902,7 +902,7 @@ "outputs": [], "source": [ "# calibrate the melt_f and annual MB \n", - "# Create a dataframe with precip factor going from 0.1 to 0.5 in 0.3 steps\n", + "# Create a dataframe with precip factor going from 0.1 to 5.0 in 0.3 steps\n", "pd_prcp_fac_sens = pd.DataFrame(index=np.arange(0.1,5.0,0.3))\n", "# Calibrate the melt factor for each of these\n", "spec_mb_prcp_fac_sens_dict = {}\n", @@ -1260,9 +1260,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Using this data-informed precipitation factor, with 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", + "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 pbservations 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 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", "\n", "![err](https://user-images.githubusercontent.com/10050469/224318400-ec1d8825-d7e7-4cdb-94f3-ebb95b8f7120.jpg)" ] @@ -1271,7 +1271,7 @@ "cell_type": "markdown", "metadata": {}, "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 imformation that we use to inform the next step.\n", + "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", "\n", @@ -1295,7 +1295,7 @@ "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 roiginal 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 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", "\n", From de11c9c45c9e88cc9b5eea00414843c84575e5e3 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 22 Aug 2024 17:22:52 +0100 Subject: [PATCH 2/3] Get ready for release (#124) --- .gitignore | 1 + _config.yml | 8 +- notebooks/10minutes/dynamical_spinup.ipynb | 2 +- notebooks/10minutes/machine_learning.ipynb | 2 +- .../10minutes/preprocessed_directories.ipynb | 13 +- notebooks/10minutes/run_with_gcm.ipynb | 2 +- .../inversion_with_frontal_ablation.ipynb | 8 +- .../kcalving_parameterization.ipynb | 454 +++++++++++++++++- .../tutorials/massbalance_global_params.ipynb | 75 +-- notebooks/welcome.ipynb | 6 +- 10 files changed, 519 insertions(+), 52 deletions(-) diff --git a/.gitignore b/.gitignore index 7cd09753..29c2789c 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ out.* RGI60-* sandbox/ outputs +ignore/ # Mac .DS_Store diff --git a/_config.yml b/_config.yml index c090a4d5..d58c73c8 100755 --- a/_config.yml +++ b/_config.yml @@ -21,19 +21,21 @@ html: use_repository_button: true use_issues_button: true use_edit_page_button: true - google_analytics_id: UA-106829797-2 extra_footer: |

These notebooks are licensed under a BSD-3-Clause license.
- © Copyright 2014-2021. + © Copyright 2014-2024.

sphinx: config: html_show_copyright: false + html_last_updated_fmt: '%b %d, %Y' nb_merge_streams: true + html_js_files: + - ['https://plausible.oggm.org/js/script.js', {'defer': 'defer', 'data-domain': 'tutorials.oggm.org'}] execute: execute_notebooks: auto # off (for tests) timeout: -1 allow_errors: true -exclude_patterns: [.virtual_documents/*,README.md,sandbox/*] +exclude_patterns: [.virtual_documents/*,README.md,sandbox/*,ignore/*] diff --git a/notebooks/10minutes/dynamical_spinup.ipynb b/notebooks/10minutes/dynamical_spinup.ipynb index c54310fe..73096043 100644 --- a/notebooks/10minutes/dynamical_spinup.ipynb +++ b/notebooks/10minutes/dynamical_spinup.ipynb @@ -324,7 +324,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.12.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, diff --git a/notebooks/10minutes/machine_learning.ipynb b/notebooks/10minutes/machine_learning.ipynb index c4a720da..6bfc07fe 100644 --- a/notebooks/10minutes/machine_learning.ipynb +++ b/notebooks/10minutes/machine_learning.ipynb @@ -1039,7 +1039,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.12.4" }, "toc": { "base_numbering": 1, diff --git a/notebooks/10minutes/preprocessed_directories.ipynb b/notebooks/10minutes/preprocessed_directories.ipynb index 37e2f43b..ef3d0370 100644 --- a/notebooks/10minutes/preprocessed_directories.ipynb +++ b/notebooks/10minutes/preprocessed_directories.ipynb @@ -180,7 +180,7 @@ }, "outputs": [], "source": [ - "rgi_ids = ['RGI60-11.01328', 'RGI60-11.00897'] " + "rgi_ids = ['RGI60-11.01328', 'RGI60-11.01450'] " ] }, { @@ -189,9 +189,10 @@ "source": [ "You can provide any number of glacier identifiers to OGGM. In this case, we chose: \n", "- `RGI60-11.01328`: [Unteraar Glacier](https://en.wikipedia.org/wiki/Unteraargletscher) in the Swiss Alps\n", - "- `RGI60-11.00897`: [Hintereisferner](http://acinn.uibk.ac.at/research/ice-and-climate/projects/hintereisferner) in the Austrian Alps.\n", + "- `RGI60-11.01450`: [Aletsch Glacier](https://en.wikipedia.org/wiki/Aletsch_Glacier) in the Swiss Alps\n", "\n", "Here is a list of other glaciers you might want to try out:\n", + "- `RGI60-11.00897`: [Hintereisferner](http://acinn.uibk.ac.at/research/ice-and-climate/projects/hintereisferner) in the Austrian Alps.\n", "- `RGI60-18.02342`: [Tasman Glacier](https://en.wikipedia.org/wiki/Tasman_Glacier) in New Zealand\n", "- `RGI60-11.00787`: [Kesselwandferner](https://de.wikipedia.org/wiki/Kesselwandferner) in the Austrian Alps\n", "- ... or any other glacier identifier! You can find other glacier identifiers by exploring the [GLIMS viewer](https://www.glims.org/maps/glims).\n", @@ -257,7 +258,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "- the keyword `from_prepro_level` indicates that we will start from [pre-processed directories](https://docs.oggm.org/en/stable/shop.html#pre-processed-directories), i.e. data that are already prepared by the OGGM team. In many cases you will want to start from pre-processed directories, in most case from level 3 or 5. For level 3 and above the model has already been calibrated, so you no longer need to do that yourself and can start rigth away with your simulation. Here we start from level 4 and add some data to the processing in order to demonstrate the OGGM workflow.\n", + "- the keyword `from_prepro_level` indicates that we will start from [pre-processed directories](https://docs.oggm.org/en/stable/shop.html#pre-processed-directories), i.e. data that are already prepared by the OGGM team. In many cases you will want to start from pre-processed directories, in most cases from level 3 or 5. For level 3 and above the model has already been calibrated, so you no longer need to do that yourself and can start rigth away with your simulation. Here we start from level 4 and add some data to the processing in order to demonstrate the OGGM workflow.\n", "- the `prepro_border` keyword indicates the number of grid points which we'd like to add to each side of the glacier for the local map: the larger the glacier will grow, the larger the border parameter should be. The available pre-processed border values are: **10, 80, 160, 240** (depending on the model set-ups there might be more or less options). These are the fixed map sizes we prepared for you - any other map size will require a full processing (see the [further DEM sources example](../tutorials/dem_sources.ipynb) for a tutorial)." ] }, @@ -418,8 +419,8 @@ "with xr.open_dataset(gdir.get_filepath('climate_historical')) as ds:\n", " ds = ds.load()\n", "# Plot the data\n", - "ds.temp.resample(time='AS').mean().plot(label=f'Annual temperature at {int(ds.ref_hgt)}m a.s.l.');\n", - "ds.temp.resample(time='AS').mean().rolling(time=31, center=True, min_periods=15).mean().plot(label='30yr average');\n", + "ds.temp.resample(time='YS').mean().plot(label=f'Annual temperature at {int(ds.ref_hgt)}m a.s.l.');\n", + "ds.temp.resample(time='YS').mean().rolling(time=31, center=True, min_periods=15).mean().plot(label='30yr average');\n", "plt.legend();" ] }, @@ -560,7 +561,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.12.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, diff --git a/notebooks/10minutes/run_with_gcm.ipynb b/notebooks/10minutes/run_with_gcm.ipynb index 5ca7867f..3f69fc8c 100644 --- a/notebooks/10minutes/run_with_gcm.ipynb +++ b/notebooks/10minutes/run_with_gcm.ipynb @@ -438,7 +438,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.12.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, diff --git a/notebooks/construction/inversion_with_frontal_ablation.ipynb b/notebooks/construction/inversion_with_frontal_ablation.ipynb index 856aa153..7bbbc39e 100644 --- a/notebooks/construction/inversion_with_frontal_ablation.ipynb +++ b/notebooks/construction/inversion_with_frontal_ablation.ipynb @@ -4,14 +4,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Ice thickness inversion with frontal ablation: a case study" + "# Ice thickness inversion with frontal ablation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "This notebook has been temporarily removed from the tutorials, as we are working on a better version! Check our progress [here]()." + "This notebook has been temporarily removed from the tutorials.\n", + "\n", + "The notebook worked with OGGM v1.5.3: " ] }, { @@ -43,7 +45,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, diff --git a/notebooks/construction/kcalving_parameterization.ipynb b/notebooks/construction/kcalving_parameterization.ipynb index 53aeb277..9d2dc540 100644 --- a/notebooks/construction/kcalving_parameterization.ipynb +++ b/notebooks/construction/kcalving_parameterization.ipynb @@ -11,7 +11,455 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This notebook has been temporarily removed from the tutorials, as we are working on a better version! Check our progress [here]()." + "
\n", + " \n", + " The calving parameterization in OGGM has been developped for versions of OGGM before 1.6 (mostly: 1.5.3). The v1.6 series brought several changes in the mass balance calibration which made a lot of the code obsolete and in need of updates. As of now (Feb 2024), calving dynamics are implemented in the dynamical model and \"work\". There is still quite some work to do make calving (including calibration) fully operational. \n", + " \n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We implement the very simple frontal ablation parameterization of [Oerlemans & Nick (2005)](https://www.cambridge.org/core/journals/annals-of-glaciology/article/minimal-model-of-a-tidewater-glacier/C6B72F547D8C44CDAAAD337E1F2FC97F) in OGGM's flowline model:\n", + "\n", + "$$F_{calving} = k d H_f w$$\n", + "\n", + "With $H_f$, $d$ and $w$ the ice thickness, the water depth and the glacier width at the calving front, $F_{calving}$ the calving flux in units m$^3$ yr$^{-1}$ and $k$ a scaling parameter that needs to be calibrated (units: yr$^{-1}$). Another way to see at this equation is to note that the frontal calving rate $\\frac{F_{Calving}}{H_f w}$ (m yr$^{-1}$) is simply the water depth scaled by the $k$ parameter.\n", + "\n", + "[Malles et al., 2023](https://doi.org/10.1017/jog.2023.19) introduces some changes to the physics and numerics of calving. We illustrate both models below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import xarray as xr\n", + "import time, os\n", + "\n", + "from oggm.core.massbalance import ScalarMassBalance\n", + "from oggm import cfg, utils, workflow, tasks, graphics\n", + "from oggm.tests.funcs import bu_tidewater_bed\n", + "from oggm.core.flowline import FluxBasedModel\n", + "\n", + "cfg.initialize(logging_level='WARNING')\n", + "cfg.PARAMS['cfl_number'] = 0.01 # less numerical instabilities\n", + "cfg.PARAMS['use_multiprocessing'] = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## FluxBasedModel \"default\" implementation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The implementation in OGGM's flowline model was relatively straightforward but required some tricks. You can have a look at the code [here](https://github.com/OGGM/oggm/blob/b5934589466bc5245303ce762572f2df8cf92369/oggm/core/flowline.py#L1322-L1379), but in short:\n", + "- the terminus thickness $H_f$ is defined as the last grid point on the calving flowline with its surface elevation **above the water level**. That is, if small chunks of ice after that are below water, their thickness is not used for the calving equation.\n", + "- the calving flux computed with the equation above is added to a \"bucket\". This bucket can be understood as \"ice that has calved but has not yet been removed from the flowline geometry\". We remove this bucket from the total flowline volume when computing model output, though.\n", + "- when there is ice below water (e.g. due to ice deformation at the front) and the bucket is large enough, remove it and subtract its volume from the bucket.\n", + "- if, after that, the bucket is larger than the total volume of the last flowline grid point above water (the calving front), remove the calving front (calve it) and subtract its volume from the bucket.\n", + "\n", + "To avoid numerical difficulties, we introduce a \"flux limiter\" at the glacier terminus. The slope between the last grid point above water (the calving front) and the next grid point (often: the sea bed) is cropped to a maximum threshold (per default: the difference between the calving front altitude and the water level) in order to limit ice deformation. See the example below for details." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Idealized experiments " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use the idealized bed profile from [Oerlemans & Nick (2005)](https://www.cambridge.org/core/journals/annals-of-glaciology/article/minimal-model-of-a-tidewater-glacier/C6B72F547D8C44CDAAAD337E1F2FC97F) and [Bassis & Ultee (2019)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JF005160). This profile has a deepening followed by a bump, allowing to study quite interesting water-depth – calving-rate feedbacks:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bu_fl = bu_tidewater_bed()[0]\n", + "\n", + "xc = bu_fl.dis_on_line * bu_fl.dx_meter / 1000\n", + "f, ax = plt.subplots(1, 1, figsize=(12, 5))\n", + "ax.plot(xc, bu_fl.bed_h, color='k')\n", + "plt.hlines(0, *xc[[0, -1]], color='C0', linestyles=':')\n", + "plt.ylim(-350, 1000); plt.ylabel('Altitude [m]'); plt.xlabel('Distance along flowline [km]');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create a simple experiment where the surface mass-balance is zero in the entire domain. Ice enters via a flux gate on the left (unit: m$^{3}$ s$^{-1}$). We can vary the amount of ice entering the domain via the flux gate. We let the glacier grow until equilibrium with three different flux values:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mb_model = ScalarMassBalance()\n", + "\n", + "to_plot = None\n", + "keys = []\n", + "for flux_gate in [0.06, 0.10, 0.16]:\n", + " model = FluxBasedModel(bu_tidewater_bed(), mb_model=mb_model,\n", + " is_tidewater=True, \n", + " flux_gate=flux_gate, # default is 0\n", + " calving_k=0.2, # default is 2.4\n", + " do_kcalving=True\n", + " )\n", + " # long enough to reach approx. equilibrium \n", + " ds = model.run_until_and_store(6000)\n", + " df_diag = model.get_diagnostics()\n", + " \n", + " if to_plot is None:\n", + " to_plot = df_diag\n", + " \n", + " key = 'Flux gate={:.02f}. Calving rate: {:.0f} m yr-1'.format(flux_gate, model.calving_rate_myr)\n", + " to_plot[key] = df_diag['surface_h']\n", + " keys.append(key)\n", + " \n", + " # Plot of volume\n", + " (ds.volume_m3 * 1e-9).plot(label=key);\n", + "plt.legend(); plt.ylabel('Volume [km$^{3}$]');\n", + "to_plot.index = xc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f, ax = plt.subplots(1, 1, figsize=(12, 5))\n", + "to_plot[keys].plot(ax=ax);\n", + "to_plot.bed_h.plot(ax=ax, color='k')\n", + "plt.hlines(0, *xc[[0, -1]], color='C0', linestyles=':')\n", + "plt.ylim(-350, 1000); plt.ylabel('Altitude [m]'); plt.xlabel('Distance along flowline [km]');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The larger the incoming ice flux, the bigger and longer the glacier. If the flux is large enough, the glacier will grow past the deepening and the bump, until the water depth and the calving rate become large enough to compensate for the total ice entering the domain." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Flux limiter " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The frontal \"free board\" (the height of the ice above water at the terminus) is quite high in the plot above. This is because we use a \"flux limiter\", effectively reducing the shear stress at the steep glacier front and avoiding high velocities, simplifying the numerics. This limiter is of course arbitrary and non-physical, but we still recommend to switch it on for your runs. Here we show the differences between two runs (with and without flux limiter):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "to_plot = None\n", + "keys = []\n", + "for limiter in [True, False]:\n", + " model = FluxBasedModel(bu_tidewater_bed(), mb_model=mb_model,\n", + " is_tidewater=True, \n", + " calving_use_limiter=limiter, # default is True\n", + " flux_gate=0.06, # default is 0\n", + " calving_k=0.2, # default is 2.4\n", + " do_kcalving=True\n", + " )\n", + " # long enough to reach approx. equilibrium \n", + " ds = model.run_until_and_store(7000)\n", + " df_diag = model.get_diagnostics()\n", + " \n", + " if to_plot is None:\n", + " to_plot = df_diag\n", + " \n", + " key = 'Flux limiter={}. Calving rate: {:.0f} m yr-1'.format(limiter, model.calving_rate_myr)\n", + " to_plot[key] = df_diag['surface_h']\n", + " keys.append(key)\n", + " \n", + " # Plot of volume\n", + " (ds.volume_m3 * 1e-9).plot(label=key);\n", + "plt.legend(); plt.ylabel('Volume [km$^{3}$]');\n", + "to_plot.index = xc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f, ax = plt.subplots(1, 1, figsize=(12, 5))\n", + "to_plot[keys].plot(ax=ax);\n", + "to_plot.bed_h.plot(ax=ax, color='k')\n", + "plt.hlines(0, *xc[[0, -1]], color='C0', linestyles=':')\n", + "plt.ylim(-350, 1000); plt.ylabel('Altitude [m]'); plt.xlabel('Distance along flowline [km]');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the rest of this notebook, we will keep the flux limiter **on**. It is also switched on per default in OGGM." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# see\n", + "cfg.PARAMS['calving_use_limiter']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Water-depth – calving-rate feedbacks" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's have some fun! We apply a periodic forcing to our glacier and study the advance and retreat of our glacier (the simulation below can take a couple minutes to run)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "years = np.arange(6001)\n", + "flux = 0.4 + 0.4 * np.sin(2 * np.pi * years / 5000)\n", + "def flux_gate(year):\n", + " return flux[int(year)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model = FluxBasedModel(bu_tidewater_bed(), mb_model=mb_model,\n", + " is_tidewater=True, \n", + " glen_a=cfg.PARAMS['glen_a']*3, # make the glacier flow faster\n", + " flux_gate=flux_gate, # default is 0\n", + " calving_k=1, # default is 2.4\n", + " do_kcalving=True\n", + " )\n", + "t0 = time.time()\n", + "ds = model.run_until_and_store(len(flux)-1)\n", + "print('Done! Time needed: {}s'.format(int(time.time()-t0)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Prepare the data for plotting\n", + "df = (ds.volume_m3 * 1e-9).to_dataframe(name='Volume [km$^3$]')[['Volume [km$^3$]']]\n", + "df['Length [m]'] = (ds['length_m'] / 1000).to_series()\n", + "df['Calving rate [m y$^{-1}$]'] = ds['calving_rate_myr'].to_series()\n", + "df['Forcing'] = flux\n", + "\n", + "# Thresholds\n", + "deep_val = 27\n", + "dfs = df.loc[(df['Length [m]'] >= deep_val) & (df.index < 5000)]\n", + "deep_t0, deep_t1 = dfs.index[0], dfs.index[-1]\n", + "dfs = df.loc[(df['Length [m]'] >= deep_val) & (df.index > 5000)]\n", + "deep_t2 = dfs.index[0]\n", + "\n", + "bump_val = 37.5\n", + "dfs = df.loc[(df['Length [m]'] >= bump_val) & (df.index < 5000)]\n", + "bump_t0, bump_t1 = dfs.index[0], dfs.index[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# The plot\n", + "f, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(9, 9), sharex=True)\n", + "ts = df['Forcing']\n", + "ts.plot(ax=ax1, color='C0');\n", + "ax1.set_ylabel(ts.name)\n", + "ts = df['Length [m]']\n", + "ts.plot(ax=ax2, color='C1');\n", + "ax2.hlines(deep_val, deep_t0, deep_t1, color='black', linestyles=':')\n", + "ax2.hlines(deep_val, deep_t2, 6000, color='black', linestyles=':')\n", + "ax2.hlines(bump_val, bump_t0, bump_t1, color='grey', linestyles='--')\n", + "ax2.annotate('Deepening', (deep_t0, deep_val-5))\n", + "ax2.annotate('Bump', (bump_t0, bump_val-5))\n", + "ax2.set_ylabel(ts.name)\n", + "# The calving rate is a bit noisy because of the bucket trick - we smooth\n", + "ts = df['Calving rate [m y$^{-1}$]'].rolling(11, center=True).max()\n", + "ts.plot(ax=ax3, color='C3')\n", + "ax3.vlines([deep_t0, deep_t1, deep_t2], ts.min(), ts.max(), color='black', linestyles=':')\n", + "ax3.vlines([bump_t0, bump_t1], ts.min(), ts.max(), color='grey', linestyles='--');\n", + "ax3.set_ylabel(ts.name); ax3.set_xlabel('Years');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our simple model reproduces the results of [Oerlemans & Nick (2005)](https://www.cambridge.org/core/journals/annals-of-glaciology/article/minimal-model-of-a-tidewater-glacier/C6B72F547D8C44CDAAAD337E1F2FC97F) and [Bassis & Ultee (2019)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JF005160) qualitatively (with the models and model parameters being different, an exact comparison is difficult). For example, here is Fig. 8 from [Bassis & Ultee (2019)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JF005160):" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Same experiments, but with the new implementation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Malles et al., 2023](https://doi.org/10.1017/jog.2023.19) uses the same equation but makes substantial changes to some of the physics. In a nutshell:\n", + "- backpressure from the water is implemented, making physics more realistic and removing the need for the flux limiter\n", + "- sliding is increased when below water\n", + "- some other logics and checks that where needed to make calving more realistic in edge cases (this is a large part of the code)\n", + "\n", + "Otherwise it works pretty much like the existing one:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from oggm.sandbox import calving_jan" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model = calving_jan.CalvingFluxBasedModelJan(bu_tidewater_bed(), mb_model=mb_model,\n", + " is_tidewater=True, \n", + " glen_a=cfg.PARAMS['glen_a']*3, # make the glacier flow faster\n", + " flux_gate=flux_gate, # default is 0\n", + " calving_k=0.6, # default is 2.4 - Note the different value for calving here to obtain similar results as above\n", + " do_kcalving=True\n", + " )\n", + "t0 = time.time()\n", + "ds_new = model.run_until_and_store(len(flux)-1)\n", + "print('Done! Time needed: {}s'.format(int(time.time()-t0)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First observation is that its slower by a factor at least two. Let's have a look at the results:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Prepare the data for plotting\n", + "df = (ds.volume_m3 * 1e-9).to_dataframe(name='Volume [km$^3$]')[['Volume [km$^3$]']]\n", + "df['Length [m]'] = (ds_new['length_m'] / 1000).to_series()\n", + "df['Calving rate [m y$^{-1}$]'] = ds_new['calving_rate_myr'].to_series()\n", + "df['Forcing'] = flux\n", + "\n", + "# Thresholds\n", + "deep_val = 27\n", + "dfs = df.loc[(df['Length [m]'] >= deep_val) & (df.index < 5000)]\n", + "deep_t0, deep_t1 = dfs.index[0], dfs.index[-1]\n", + "dfs = df.loc[(df['Length [m]'] >= deep_val) & (df.index > 5000)]\n", + "deep_t2 = dfs.index[0]\n", + "\n", + "bump_val = 37.5\n", + "dfs = df.loc[(df['Length [m]'] >= bump_val) & (df.index < 5000)]\n", + "bump_t0, bump_t1 = dfs.index[0], dfs.index[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# The plot\n", + "f, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(9, 9), sharex=True)\n", + "ts = df['Forcing']\n", + "ts.plot(ax=ax1, color='C0');\n", + "ax1.set_ylabel(ts.name)\n", + "ts = df['Length [m]']\n", + "ts.plot(ax=ax2, color='C1');\n", + "ax2.hlines(deep_val, deep_t0, deep_t1, color='black', linestyles=':')\n", + "ax2.hlines(deep_val, deep_t2, 6000, color='black', linestyles=':')\n", + "ax2.hlines(bump_val, bump_t0, bump_t1, color='grey', linestyles='--')\n", + "ax2.annotate('Deepening', (deep_t0, deep_val-5))\n", + "ax2.annotate('Bump', (bump_t0, bump_val-5))\n", + "ax2.set_ylabel(ts.name)\n", + "# The calving rate is a bit noisy because of the bucket trick - we smooth\n", + "ts = df['Calving rate [m y$^{-1}$]'].rolling(11, center=True).max()\n", + "ts.plot(ax=ax3, color='C3')\n", + "ax3.vlines([deep_t0, deep_t1, deep_t2], ts.min(), ts.max(), color='black', linestyles=':')\n", + "ax3.vlines([bump_t0, bump_t1], ts.min(), ts.max(), color='grey', linestyles='--');\n", + "ax3.set_ylabel(ts.name); ax3.set_xlabel('Years'); ax3.set_ylim(0, 150);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds.length_m.plot();\n", + "ds_new.length_m.plot();" ] }, { @@ -21,7 +469,7 @@ "## What's next?\n", "\n", "- return to the [OGGM documentation](https://docs.oggm.org)\n", - "- back to the [table of contents](../welcome.ipynb)" + "- back to the [table of contents](welcome.ipynb)" ] } ], @@ -42,7 +490,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.4" }, "toc": { "base_numbering": 1, diff --git a/notebooks/tutorials/massbalance_global_params.ipynb b/notebooks/tutorials/massbalance_global_params.ipynb index 62ed6fe4..3553f361 100644 --- a/notebooks/tutorials/massbalance_global_params.ipynb +++ b/notebooks/tutorials/massbalance_global_params.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "2ff429f6-c55f-4e1a-a866-a62b00cdfabd", + "id": "0", "metadata": {}, "source": [ "# Global distribution of the mass-balance model parameters" @@ -10,7 +10,7 @@ }, { "cell_type": "markdown", - "id": "8fa17770-c3f5-47aa-8fa9-083e49c9c13d", + "id": "1", "metadata": {}, "source": [ "This notebook is a follow-up on the very important overview of the [mass-balance calibration procedure in v1.6](massbalance_calibration.ipynb).\n", @@ -24,7 +24,7 @@ { "cell_type": "code", "execution_count": null, - "id": "74bbe161-f091-48cf-91c9-9ad183fc35ab", + "id": "2", "metadata": {}, "outputs": [], "source": [ @@ -37,7 +37,7 @@ }, { "cell_type": "markdown", - "id": "aa5a7511-cd6d-48b8-aaea-184ee6bf050e", + "id": "3", "metadata": {}, "source": [ "Let's start by reading in the \"glacier statistics\" files, a bunch of statistics written by OGGM at the end of the preprocessing:" @@ -46,7 +46,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8235658d-c668-46fb-8af2-df63df5dd464", + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -68,7 +68,7 @@ }, { "cell_type": "markdown", - "id": "bbb32744-cb85-4eee-852e-f2283a0e2e93", + "id": "5", "metadata": {}, "source": [ "## Global statistics of OGGM's 1.6.1 \"informed three steps\" method" @@ -76,7 +76,7 @@ }, { "cell_type": "markdown", - "id": "2a3c9ab8-cceb-40fc-884a-fff64eb5fcc8", + "id": "6", "metadata": {}, "source": [ "As explained in the [mass-balance calibration procedure in v1.6](massbalance_calibration.ipynb) notebook, the \"informed three steps\" method provides first guesses for the precipitation factor and the temperature bias. We then calibrate each glacier in three steps - let's check the number of glaciers calibrated this way:" @@ -84,7 +84,7 @@ }, { "cell_type": "markdown", - "id": "0193adb7-4543-4ebe-9597-c53e69435b27", + "id": "7", "metadata": {}, "source": [ "Step 0: use the first guess `melt_f` = 5, `prcp_fac` = data-informed from winter precipitation, `temp_bias` = data-informed from the global calibration with fixed parameters (see [mass-balance calibration procedure in v1.6](massbalance_calibration.ipynb) for details). \n", @@ -97,7 +97,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5a09dc13-661a-4a48-8bcb-b5a89cf4636e", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -109,7 +109,7 @@ }, { "cell_type": "markdown", - "id": "6136d2c6-b8f6-424a-a7a0-85f5323260cd", + "id": "9", "metadata": {}, "source": [ "Step 2: if Step 1 did not work, we allow `melt_f` to vary between a predefined range (1.5 - 17) while fixing `temp_bias` and `prcp_fac` again.\n", @@ -122,7 +122,7 @@ { "cell_type": "code", "execution_count": null, - "id": "611601be-b346-4b0f-b727-a60af9fc159c", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -140,7 +140,7 @@ }, { "cell_type": "markdown", - "id": "9d29e23b-7c56-4e56-af1a-e512835b1ee6", + "id": "11", "metadata": {}, "source": [ "## Global parameter distributions" @@ -148,7 +148,7 @@ }, { "cell_type": "markdown", - "id": "d999d5c2-980e-48b5-8485-950862167066", + "id": "12", "metadata": {}, "source": [ "### Melt factor" @@ -157,7 +157,7 @@ { "cell_type": "code", "execution_count": null, - "id": "648943c4-5630-41db-85d4-d4dd84f26c52", + "id": "13", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +178,7 @@ }, { "cell_type": "markdown", - "id": "ff4a7a01-ef37-48ca-9299-94b1c0ed0d11", + "id": "14", "metadata": {}, "source": [ "### Precip factor" @@ -187,7 +187,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18195fcc-046c-4c60-93b9-94428d4ac8c1", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -206,9 +206,20 @@ "ax2.set_ylabel('Frequency (log scale)');" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"The precipitation factor median is {df_params['prcp_fac'].median():.1f}.\")\n", + "print(f\"The 5% percentile is {df_params['prcp_fac'].quantile(0.05):.1f} and the 95% percentile is {df_params['prcp_fac'].quantile(0.95):.1f}\")" + ] + }, { "cell_type": "markdown", - "id": "ad283a65-06ba-49e9-8bf3-d84038a5e2b0", + "id": "17", "metadata": {}, "source": [ "### Temperature bias" @@ -217,7 +228,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8feca9db-b04a-4be1-834a-3c8bd3110e73", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -238,7 +249,7 @@ }, { "cell_type": "markdown", - "id": "6b8e0d90-7e19-4285-a5be-cd70ea676b80", + "id": "19", "metadata": {}, "source": [ "### Take home" @@ -246,7 +257,7 @@ }, { "cell_type": "markdown", - "id": "b1628f65-f155-4911-86d7-4044e771ee46", + "id": "20", "metadata": {}, "source": [ "- a substantial (33%) part of all glaciers are attributed the default melt factor of 5 after the first guesses in climate data bias correction. In other words, this means that we are substantially correcting the climate forcing to \"match\" the presence of a glacier. Other calibration methods are using similar techniques (they differ in the details and the allowed range of parameter values)\n", @@ -257,7 +268,7 @@ }, { "cell_type": "markdown", - "id": "1bc9f8de-7583-4d6a-b2af-49eb0e56987a", + "id": "21", "metadata": {}, "source": [ "## Influence of dynamical spinup " @@ -265,7 +276,7 @@ }, { "cell_type": "markdown", - "id": "4d8e2e86-c8d4-4b24-b831-abe263a32bdf", + "id": "22", "metadata": {}, "source": [ "The dynamical spinup procedure (explained in this [10 minutes tutorial](../10minutes/dynamical_spinup.ipynb) and in more detail in [this tutorial](dynamical_spinup.ipynb)) starts from the parameters calibrated above with a *static* geometry and calibrate the melt factor again using an iterative procedure, making sure that the parameters and the past evolution of the glacier are consistent with the past evolution of the glacier. In doing so, it achieves two things:\n", @@ -278,7 +289,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f4143761-4d76-4e19-b1b5-7f4b2001951c", + "id": "23", "metadata": {}, "outputs": [], "source": [ @@ -296,7 +307,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e349ab56-fc49-4425-b3ff-b35420a90342", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -305,7 +316,7 @@ }, { "cell_type": "markdown", - "id": "d2afbe29-eb76-4908-8ead-42b25ac25f19", + "id": "25", "metadata": {}, "source": [ "First of all, let's see how many glaciers have had their melt factor changed as a result of the dynamical calibration (i.e. dynamical calibration was succesful):" @@ -314,7 +325,7 @@ { "cell_type": "code", "execution_count": null, - "id": "860c82b1-4812-4079-a4bc-50a27cffa6ff", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -326,7 +337,7 @@ }, { "cell_type": "markdown", - "id": "86ea2057-d2b2-48e4-aaf5-d3736dc39525", + "id": "27", "metadata": {}, "source": [ "Let's plot the change in distribution of the parameters:" @@ -335,7 +346,7 @@ { "cell_type": "code", "execution_count": null, - "id": "855acbd5-e3fe-4be0-98aa-006a4cbb6720", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -356,7 +367,7 @@ }, { "cell_type": "markdown", - "id": "36c0d60d-c8ca-4e6c-9fe2-bd175e9fd949", + "id": "29", "metadata": {}, "source": [ "In which direction is the parameter changed?" @@ -365,7 +376,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9031f258-d474-492e-a1d4-2a4b6f966df6", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -385,7 +396,7 @@ }, { "cell_type": "markdown", - "id": "b7465f5b-bc52-416d-8db4-e168a6ff89be", + "id": "31", "metadata": {}, "source": [ "### Take home\n", @@ -395,7 +406,7 @@ }, { "cell_type": "markdown", - "id": "4d2bfc42-e753-4745-b9cb-fc75d75cae41", + "id": "32", "metadata": { "jp-MarkdownHeadingCollapsed": true }, diff --git a/notebooks/welcome.ipynb b/notebooks/welcome.ipynb index 1fb6b221..40d45441 100644 --- a/notebooks/welcome.ipynb +++ b/notebooks/welcome.ipynb @@ -32,7 +32,6 @@ "- [deal_with_errors](tutorials/deal_with_errors.ipynb): dealing with errors after a run\n", "- [elevation_bands_vs_centerlines](tutorials/elevation_bands_vs_centerlines.ipynb): differences between \"elevation band\" and \"centerline\" flowlines\n", "- [full_prepro_workflow](tutorials/full_prepro_workflow.ipynb): what's in your preprocessed directories? A full OGGM workflow, step by step\n", - "- [rebuilding_the_prepro_gdirs.ipynb](tutorials/rebuilding_the_prepro_gdirs.ipynb): Step-by-Step guide to rebuilding the preprocessed directories from scratch\n", "\n", "## Mass balance\n", "\n", @@ -60,6 +59,9 @@ "- [inversion](tutorials/inversion.ipynb): run the OGGM ice thickness inversion model with various ice parameters\n", "- [observed_thickness_with_dynamic_spinup](tutorials/observed_thickness_with_dynamic_spinup.ipynb): how to create OGGM flowlines from thickness observations and dynamically initialise the model\n", "\n", + "## Calving\n", + "\n", + "\n", "## OGGM shop and additional data\n", "\n", "- [oggm_shop](tutorials/oggm_shop.ipynb): OGGM-Shop and Glacier Directories in OGGM\n", @@ -129,7 +131,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.12.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, From 6ce87a90de2defda52e35e3a670988d498618192 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 22 Aug 2024 19:03:56 +0100 Subject: [PATCH 3/3] Add key --- .github/workflows/build-pages.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/build-pages.yml b/.github/workflows/build-pages.yml index 61c23d9e..ddd82949 100644 --- a/.github/workflows/build-pages.yml +++ b/.github/workflows/build-pages.yml @@ -36,6 +36,8 @@ jobs: ${PIP} install -r requirements.txt ${PIP} uninstall -y progressbar2 - name: Build Book + env: + STATIC_MAP_API_KEY: ${{ secrets.STATIC_MAP_API_KEY }} run: | jupyter-book build . - name: Upload Build Artifacts