{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Inputs for the Probable Maximum flood (PMF)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Probable Maximum Precipitation (PMP)\n", "\n", "Probable Maximum Precipitation (PMP) is the theoretical maximum amount of precipitation that could occur at a specific location within a given period of time, considering the most extreme meteorological conditions. PMP is a critical parameter in hydrology, especially for the design of infrastructure such as dams, reservoirs, and drainage systems.\n", "\n", "There are several methods for calculating PMP, each varying in complexity and the type of data used. The method currently implemented in `xHydro` is based on the approach outlined by [Clavet-Gaumont et al. (2017)](https://doi.org/10.1016/j.ejrh.2017.07.003). This method involves maximizing the precipitable water over a given location, which refers to the total water vapor in the atmosphere that could potentially be converted into precipitation under ideal conditions. By maximizing this value, the method estimates the maximum precipitation that could theoretically occur at the location.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "import hvplot.xarray\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pooch\n", "import xarray as xr\n", "import xclim\n", "\n", "import xhydro as xh\n", "from xhydro.testing.helpers import deveraux" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Acquiring data\n", "\n", "The acquisition of climatological data is outside the scope of `xHydro`. However, some examples of how to obtain and handle such data are provided in the [GIS operations](../gis.ipynb) and [Use Case Example](../use_case.ipynb) notebooks. For this notebook, we will use a test dataset consisting of 2 years and 3x3 grid cells from CanESM5 climate model data. In a real application, it would be preferable to have as many years of data as possible.\n", "\n", "To perform the analysis, certain climatological variables are required.\n", "\n", "- **Daily Timestep Variables**:\n", " - `pr` → Precipitation flux\n", " - `snw` → Snow water equivalent\n", " - `hus` → Specific humidity for multiple pressure levels\n", " - `zg` → Geopotential height for multiple pressure levels\n", "\n", "- **Fixed Field Variables**:\n", " - `orog` → Surface altitude\n", "\n", "In cold regions, it may be necessary to split total precipitation into rainfall and snowfall components. Many climate models already provide this data separately. However, if this data is not directly available, libraries such as `xclim` can approximate the split using precipitation and temperature data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "import xhydro as xh\n", "\n", "path_day_zip = deveraux().fetch(\n", " \"pmp/CMIP.CCCma.CanESM5.historical.r1i1p1f1.day.gn.zarr.zip\",\n", " pooch.Unzip(),\n", ")\n", "ds_day = xr.open_zarr(Path(path_day_zip[0]).parents[0])\n", "\n", "path_fx_zip = deveraux().fetch(\n", " \"pmp/CMIP.CCCma.CanESM5.historical.r1i1p1f1.fx.gn.zarr.zip\",\n", " pooch.Unzip(),\n", ")\n", "ds_fx = xr.open_zarr(Path(path_fx_zip[0]).parents[0])\n", "\n", "# There are a few issues with attributes in this dataset that we need to address\n", "ds_day[\"pr\"].attrs = {\"units\": \"mm\", \"long_name\": \"precipitation\"}\n", "ds_day[\"prsn\"].attrs = {\"units\": \"mm\", \"long_name\": \"snowfall\"}\n", "ds_day[\"rf\"].attrs = {\"units\": \"mm\", \"long_name\": \"rainfall\"}\n", "\n", "# Combine both datasets\n", "ds = ds_day.convert_calendar(\"standard\")\n", "ds[\"orog\"] = ds_fx[\"orog\"]\n", "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing the PMP\n", "\n", "The method outlined by [Clavet-Gaumont et al. (2017)](https://doi.org/10.1016/j.ejrh.2017.07.003) follows these steps:\n", "\n", "1. **Identification of Major Precipitation Events**: \n", " The first step involves identifying the major precipitation events that will be maximized. This is done by filtering events based on a specified threshold.\n", "\n", "2. **Computation of Monthly 100-Year Precipitable Water**: \n", " The next step involves calculating the 100-year precipitable water on a monthly basis using the Generalized Extreme Value (GEV) distribution, with a maximum cap of 20% greater than the largest observed value.\n", "\n", "3. **Maximization of Precipitation During Events**: \n", " In this step, the precipitation events are maximized based on the ratio between the 100-year monthly precipitable water and the precipitable water during the major precipitation events. In snow-free regions, this is the final result.\n", "\n", "4. **Seasonal Separation in Cold Regions**: \n", " In cold regions, the results are separated into seasons (e.g., spring, summer) to account for snow during the computation of Probable Maximum Floods (PMF).\n", "\n", "This method provides a comprehensive approach for estimating the PMP, taking into account both temperature and precipitation variations across different regions and seasons.\n", "\n", "\n", "### Major precipitation events\n", "\n", "The first step in calculating the Probable Maximum Precipitation (PMP) involves filtering the precipitation data to retain only the events that exceed a certain threshold. These major precipitation events will be maximized in subsequent steps. The function `xh.indicators.pmp.major_precipitation_events` can be used for this purpose. It also provides the option to sum precipitation over a specified number of days, which can help aggregate storm events. For 2D data, such as in this example, each grid point is treated independently.\n", "\n", "In this example, we will filter out the 10% most intense storms to avoid overemphasizing smaller precipitation events during the maximization process. Additionally, we will focus on rainfall (`rf`) rather than total precipitation (`pr`) to exclude snowstorms and ensure that we are only considering liquid precipitation events.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(xh.indicators.pmp.major_precipitation_events)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "precipitation_events = xh.indicators.pmp.major_precipitation_events(\n", " ds.rf, windows=[1], quantile=0.9, min_prec=\"0.1 mm\"\n", ")\n", "\n", "ds.rf.isel(x=1, y=1).hvplot() * precipitation_events.isel(\n", " x=1, y=1, window=0\n", ").hvplot.scatter(color=\"red\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Daily precipitable water\n", "\n", "
WARNING\n", " \n", "This step should be avoided if possible, as it involves approximating precipitable water from the integral of specific humidity and will be highly sensitive to the number of pressure levels used. If available, users are strongly encouraged to use a variable or combination of variables that directly represent precipitable water.\n", "\n", "
\n", "\n", "Precipitable water can be estimated using `xhydro.indicators.pmp.precipitable_water` by integrating the vertical column of humidity. This process requires specific humidity, geopotential height, and elevation data. The resulting value represents the total amount of water vapor that could potentially be precipitated from the atmosphere under ideal conditions.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(xh.indicators.pmp.precipitable_water)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pw = xh.indicators.pmp.precipitable_water(\n", " hus=ds.hus,\n", " zg=ds.zg,\n", " orog=ds.orog,\n", " windows=[1],\n", " add_pre_lay=False,\n", ")\n", "\n", "pw.isel(x=1, y=1, window=0).hvplot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Monthly 100-year precipitable water\n", "\n", "According to Clavet-Gaumont et al. (2017), a monthly 100-year precipitable water must be computed using the Generalized Extreme Value (GEV) distribution. The value should be limited to a maximum of 20% greater than the largest observed precipitable water value for a given month. This approach ensures that the estimated 100-year event is realistic and constrained by observed data.\n", "\n", "To compute this, you can use the `xh.indicators.pmp.precipitable_water_100y` function. If using `rebuild_time`, the output will have the same time axis as the original data.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(xh.indicators.pmp.precipitable_water_100y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pw100 = xh.indicators.pmp.precipitable_water_100y(\n", " pw.sel(window=1).chunk(dict(time=-1)),\n", " dist=\"genextreme\",\n", " method=\"ML\",\n", " mf=0.2,\n", " rebuild_time=True,\n", ").compute()\n", "\n", "pw.isel(x=1, y=1, window=0).hvplot() * pw100.isel(x=1, y=1).hvplot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.array([2, 5, 6, 7, 7])\n", "y = [5, 7]\n", "\n", "np.isin(x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Maximized precipitation\n", "\n", "
INFO\n", " \n", "This step follows the methodology described in Clavet-Gaumont et al., 2017. It is referred to as \"Maximizing precipitation\", however, it effectively applies a ratio based on the monthly 100-year precipitable water. If a historical event surpassed this value—such as the case observed for January 2011—the result may actually lower the precipitation, rather than increasing it.\n", "\n", "
\n", "\n", "With the information gathered so far, we can now proceed to maximize the precipitation events. Although `xHydro` does not provide an explicit function for this step, it can be accomplished by following these steps:\n", "\n", "1. **Compute the Ratio**: First, calculate the ratio between the 100-year monthly precipitable water and the precipitable water during the major precipitation events.\n", " \n", "2. **Apply the Ratio**: Next, apply this ratio to the precipitation values themselves to maximize the precipitation events accordingly.\n", "\n", "This process effectively scales the precipitation events based on the 100-year precipitable water, giving an estimate of the maximum possible rainfall.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Precipitable water on the day of the major precipitation events.\n", "pw_events = pw.where(precipitation_events > 0)\n", "ratio = pw100 / pw_events\n", "\n", "# Apply the ratio onto precipitation itself\n", "precipitation_max = ratio * precipitation_events\n", "precipitation_max.name = \"maximized_precipitation\"\n", "\n", "ds.rf.isel(x=1, y=1).hvplot() * precipitation_max.isel(\n", " x=1, y=1, window=0\n", ").hvplot.scatter(color=\"red\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Seasonal Mask\n", "\n", "In cold regions, computing Probable Maximum Floods (PMFs) often involves scenarios that combine both rainfall and snowpack. Therefore, PMP values may need to be separated into two categories: rain-on-snow (i.e., \"spring\") and snow-free rainfall (i.e., \"summer\").\n", "\n", "This can be computed easily using `xhydro.indicators.pmp.compute_spring_and_summer_mask`, which defines the start and end dates of spring, summer, and winter based on the presence of snow on the ground, with the following criteria:\n", "\n", "1. **Winter**: \n", " - Winter start: The first day after which there are at least 14 consecutive days with snow on the ground. \n", " - Winter end: The last day with snow on the ground, followed by at least 45 consecutive snow-free days.\n", "\n", "2. **Spring**: \n", " - Spring start: 60 days before the end of winter.\n", " - Spring end: 30 days after the end of winter.\n", "\n", "3. **Summer**: \n", " - The summer period is defined as the time between winters. This period is not influenced by whether it falls in the traditional summer or fall seasons, but rather simply marks the interval between snow seasons.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(xh.indicators.pmp.compute_spring_and_summer_mask)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mask = xh.indicators.pmp.compute_spring_and_summer_mask(\n", " ds.snw,\n", " thresh=\"1 cm\",\n", " window_wint_end=14, # Since the dataset used does not have a lot of snow, we need to be more lenient\n", " freq=\"YS-SEP\",\n", ")\n", "\n", "mask" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xclim.core.units.convert_units_to(\n", " ds.isel(x=1, y=1).snw, \"cm\", context=\"hydro\"\n", ").hvplot() * (mask.mask_spring.isel(x=1, y=1) * 10).hvplot() * (\n", " mask.mask_summer.isel(x=1, y=1) * 8\n", ").hvplot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Final PMP\n", "\n", "The final PMP is obtained by finding the maximum value over the `time` dimension. In our case, since we computed a season mask, we can further refine the results into a spring and summer PMP." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pmp_spring = (precipitation_max * mask.mask_spring).max(\"time\").compute()\n", "pmp_summer = (precipitation_max * mask.mask_summer).max(\"time\").compute()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.subplots(1, 2, figsize=[12, 5])\n", "\n", "ax = plt.subplot(1, 2, 1)\n", "pmp_spring.sel(window=1).plot(vmin=30, vmax=100)\n", "plt.title(\"Spring PMP\")\n", "\n", "ax = plt.subplot(1, 2, 2)\n", "pmp_summer.sel(window=1).plot(vmin=30, vmax=100)\n", "plt.title(\"Summer PMP\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PMPs with aggregated storm configurations\n", "\n", "In some cases, it may be preferable to avoid processing each grid cell independently. Instead, storms can be aggregated using various configurations to provide a more regionally representative estimate. These configurations allow for the spatial averaging of storm events, which can help reduce variability across grid cells and yield more reliable results.\n", "\n", "Different aggregation configurations are discussed in Clavet-Gaumont et al. (2017) and have been implemented in `xHydro` under the function `xhydro.indicators.pmp.spatial_average_storm_configurations`.\n", "\n", "Note that precipitable water must first be calculated in a distributed manner and then spatially averaged to obtain the aggregated precipitable water.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(xh.indicators.pmp.spatial_average_storm_configurations)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds_agg = []\n", "for variable in [\"rf\", \"pw\", \"snw\"]:\n", " if variable == \"pw\":\n", " ds_agg.extend(\n", " [xh.indicators.pmp.spatial_average_storm_configurations(pw, radius=3)]\n", " )\n", " else:\n", " ds_agg.extend(\n", " [\n", " xh.indicators.pmp.spatial_average_storm_configurations(\n", " ds[variable], radius=3\n", " )\n", " ]\n", " )\n", "ds_agg = xr.merge(ds_agg).chunk(dict(time=-1))\n", "\n", "# The aggreagtion creates NaN values for snow, so we'll restrict the domain\n", "ds_agg = ds_agg.isel(y=slice(0, -1), x=slice(0, -1))\n", "\n", "ds_agg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After applying storm aggregation, the subsequent steps remain the same as before, following the standard PMP calculation process outlined earlier." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pe_agg = xh.indicators.pmp.major_precipitation_events(\n", " ds_agg.rf, windows=[1], quantile=0.9, min_prec=\"0.1 mm\"\n", ")\n", "\n", "pw100_agg = xh.indicators.pmp.precipitable_water_100y(\n", " ds_agg.sel(window=1).precipitable_water, dist=\"genextreme\", method=\"ML\", mf=0.2\n", ")\n", "\n", "# Maximization\n", "pw_events_agg = ds_agg.precipitable_water.where(pe_agg > 0)\n", "r_agg = pw100_agg / pw_events_agg\n", "\n", "pmax_agg = r_agg * pe_agg\n", "\n", "# Season mask\n", "mask_agg = xh.indicators.pmp.compute_spring_and_summer_mask(\n", " ds_agg.snw,\n", " thresh=\"1 cm\",\n", " window_wint_start=14,\n", " window_wint_end=14,\n", " spr_start=60,\n", " spr_end=30,\n", " freq=\"YS-SEP\",\n", ")\n", "\n", "pmp_spring_agg = pmax_agg * mask_agg.mask_spring\n", "pmp_summer_agg = pmax_agg * mask_agg.mask_summer\n", "\n", "pmp_summer_agg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Previously, the final PMP for each season was obtained by taking the maximum value over the `time` dimension. In this updated approach, we can now take the maximum across both the `time` and `conf` dimensions, using our multiple storm configurations.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Final results\n", "pmp_spring_agg = pmp_spring_agg.max(dim=[\"time\", \"conf\"])\n", "pmp_summer_agg = pmp_summer_agg.max(dim=[\"time\", \"conf\"])\n", "\n", "pmp_summer_agg" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Compare results for the central grid cell\n", "print(\n", " f\"Grid-cell summer PMP: {np.round(pmp_summer.isel(x=1, y=1, window=0).values, 1)} mm\"\n", ")\n", "print(\n", " f\"Aggregated summer PMP: {np.round(pmp_summer_agg.isel(x=1, y=1, window=0).values, 1)} mm\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Probable Maximum Snow Accumulation (PMSA)\n", "\n", "The PMSA represents the theoretical maximum snow water equivalent (SWE) that could accumulate in a year under the most extreme meteorological conditions. The method currently implemented in `xHydro` follows the approach introduced by [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031). In this study, the authors developed a PMSA estimation technique using simulated data from regional climate models to maximize the precipitable water leading to snowfall.\n", "\n", "\n", "The method outlined by [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) follows these steps:\n", "\n", "1. **Identification of the Precipitable Water leading to Snowfall**: \n", " Precipitable water data are filtered using specified thresholds for rainfall and snowfall to isolate values associated with snowfall events.\n", "\n", "2. **Computation of Monthly 100-Year Precipitable Water leading to Snowfall**: \n", " Monthly 100-year precipitable water values associated with snowfall are estimated using the Gamma distribution.\n", "\n", "3. **Maximization of Snowfall Events**: \n", " Snowfall events are maximized based on the ratio between the 100-year monthly precipitable water leading to snowfall and the precipitable water observed during snowfall events.\n", "\n", "4. **PMSA estimation**: \n", " For each winter, both maximized and non-maximized snowfall events are aggregated. The PMSA is then defined as the highest total among all winters.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Identification of the Precipitable Water leading to Snowfall\n", "\n", "[Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) proposed that only precipitable water values from time steps with at least 0.25 mm/6h of solid precipitation (expressed as water equivalent) and less than 0.1 mm/6h of rain should be considered. Since in this example we are working with daily time steps, these thresholds were adjusted to 1 mm/day and 0.4 mm/day, respectively." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prsn_threshold = \"1 mm\"\n", "prra_threshold = \"0.4 mm\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some cases, both snowfall and rainfall occur during the same time step. However, it is typically not possible to distinguish between them. To avoid including precipitable water values associated with rainfall events, the authors proposed the following three methods:\n", "\n", " - **m1**: Only snowfall events that are not accompanied by rain, are considered for maximization. For the precipitable water corresponding to a certain event, we consider only time steps that show at least the minimum amount of snowfall and less than the minimum amount of rain.\n", " - **m2**: Snowfall events that exceed the established threshold for maximization (i.e., minimum amount of snowfall) are maximized, regardless whether they are accompanied by rain or not. For the precipitable water corresponding to a certain event, we consider only time steps that show at least the minimum amount snowfall.\n", " - **m3**: Same events as those selected for M2, but if additionally more than the minimum amount of rain occurs, the precipitable water is multiplied by the ratio of snowfall [mm water equivalent] over total precipitation. The multiplication provides a means to estimate the amount of precipitable water that leads to snowfall.\n", "\n", "We start by calculating the precipitation events and their rain and snow components for a given accumulation. \n", "\n", "
INFO\n", " \n", "To keep the example simple, the duration of the precipitation event is set to one time step (windows=[1]). When computing the PMSA, users should pay close attention to the window parameter as it can significantly influences the PMSA results. [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) suggests to use an algorithm that actually determines the real duration of every simulated snowfall rather than using predefined durations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "precipitation_events = xh.indicators.pmp.major_precipitation_events(ds.pr, windows=[1], min_prec=\"0.1 mm\")\n", "prra_events = xh.indicators.pmp.major_precipitation_events(ds.rf, windows=[1], min_prec=\"0.1 mm\")\n", "prsn_events = xh.indicators.pmp.major_precipitation_events(ds.prsn, windows=[1], min_prec=\"0.1 mm\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pw_snowfall_m1 = xh.indicators.pmp.pw_snowfall(\n", " pw,\n", " method=\"m1\",\n", " prsn_events=prsn_events,\n", " prsn_threshold=prsn_threshold,\n", " prra_events=prra_events,\n", " prra_threshold=prra_threshold,\n", ")\n", "pw_snowfall_m2 = xh.indicators.pmp.pw_snowfall(\n", " pw, method=\"m2\", prsn_events=prsn_events, prsn_threshold=prsn_threshold\n", ")\n", "pw_snowfall_m3 = xh.indicators.pmp.pw_snowfall(\n", " pw,\n", " method=\"m3\",\n", " prsn_events=prsn_events,\n", " prsn_threshold=prsn_threshold,\n", " prra_events=prra_events,\n", " prra_threshold=prra_threshold,\n", " pr_events=precipitation_events,\n", ")\n", "\n", "pw.isel(x=1, y=1, window=0).hvplot(label=\"pw\") * pw_snowfall_m1.isel(\n", " x=1, y=1, window=0\n", ").hvplot.scatter(color=\"purple\", marker=\"o\", size=75, label=\"M1\") * pw_snowfall_m2.isel(\n", " x=1, y=1, window=0\n", ").hvplot.scatter(\n", " color=\"green\", marker=\"v\", size=50, label=\"M2\"\n", ") * pw_snowfall_m3.isel(\n", " x=1, y=1, window=0\n", ").hvplot.scatter(\n", " color=\"red\", marker=\"*\", label=\"M3\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As suggested by [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031), in this example we used the methodology M3." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computation of Monthly 100-Year Precipitable Water leading to Snowfall (pw100)\n", "\n", "[Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) recommend fitting a non-stationary Gamma distribution. Currently, only the stationary Gamma distribution is implemented, but a non-stationary version is planned for a future release. The GEV distribution can also be used for extrapolation purposes. A non-stationnary version of it is available in the xhydro package.\n", "\n", "[Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) also suggests using the maximum observed precipitable water as pw100 when there are fewer than 20 data points in a month (i.e., `n=20`). Furthermore, in line with [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) methodology, pw100 is not limited (i.e., `mf=None`). However, this assumption should be reviewed for each specific case. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pw100_snow_events_m3 = xh.indicators.pmp.precipitable_water_100y(\n", " pw_snowfall_m3.sel(window=1).chunk(dict(time=-1)),\n", " dist=\"gamma\",\n", " method=\"ML\",\n", " mf=None,\n", " n=20,\n", " rebuild_time=True,\n", ").compute()\n", "\n", "pw.isel(x=0, y=1, window=0).hvplot(label=\"pw\") * pw100_snow_events_m3.isel(\n", " x=0, y=1\n", ").hvplot(label=\"pw₁₀₀\") * pw_snowfall_m3.isel(x=0, y=1, window=0).hvplot.scatter(\n", " color=\"red\", marker=\"*\", label=\"m3\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Maximization of Snowfall Events\n", "\n", "With the information gathered so far, we can now proceed to maximize the snowfall events. Although `xHydro` does not provide an explicit function for this step, it can be accomplished by following these steps:\n", "\n", "1. **Compute the maximitation ratio (r)**: Calculate the ratio between the 100-year monthly precipitable water leading to snowfall and the precipitable water during snowfall events.\n", "\n", "2. **Limit the ratio (r)**: Limit the ratio to 1≤r≤2.5 as suggested by [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031).\n", " \n", "3. **Apply the ratio**:\n", "\n", " [Klein et al. (2016)](https://doi.org/10.1016/j.jhydrol.2016.03.031) suggest maximizing only those snowfall events that exceed a given threshold for a specified duration. Therefore, only events that meet this criterion are multiplied by the corresponding maximization ratio (r). If the criterion is not met, the snowfall event is not maximized.\n", "\n", " The authors used different thresholds depending on the duration of the event:\n", "\n", " - 6-hour duration: 3, 4, 5, 6 mm\n", " - 12-hour duration: 4, 5, 6, 7 mm\n", " - 24-hour duration: 5, 6, 7, 8 mm\n", "\n", " In this example, a 6-mm threshold was selected for the 1-day duration.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "limit = 2.5\n", "maximization_threshold = \"6 mm\"\n", "maximization_threshold_converted = xclim.core.units.convert_units_to(\n", " maximization_threshold, prsn_events\n", ") # Convert thresholds to match the data's units\n", "\n", "r = pw100_snow_events_m3 / pw_snowfall_m3\n", "r_limited = xr.where(r > limit, limit, r)\n", "r_limited = xr.where(r_limited < 1, 1, r_limited)\n", "\n", "max_snow_events = xr.where(\n", " prsn_events >= maximization_threshold_converted, r_limited * prsn_events, np.nan\n", ")\n", "max_snow_events.name = \"max_snow_events\"\n", "non_max_snow_events = prsn_events.where(prsn_events < maximization_threshold_converted)\n", "non_max_snow_events.name = \"non_max_snow_events\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximized snow accumulation (MSA) is the sum of the maximized and non-maximized events throughout a winter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "snow_sum = xr.concat([non_max_snow_events, max_snow_events], dim=\"variable\").sum(\n", " \"variable\"\n", ")\n", "\n", "snow_sum.isel(x=1, y=1, window=0).hvplot(\n", " label=\"Maximized snowfall\", color=\"red\"\n", ") * max_snow_events.isel(x=1, y=1, window=0).hvplot.scatter(\n", " color=\"red\", label=\"Maximized events\"\n", ") * prsn_events.isel(\n", " x=1, y=1, window=0\n", ").hvplot(\n", " label=\"Snowfall\", color=\"#30a2da\", line_dash=\"dashed\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PMSA estimation\n", "\n", "The PMSA is the highest total among all winters defined." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pmsa = snow_sum.resample(time=\"YS-JUL\").sum(dim=\"time\").max(dim=\"time\")\n", "pmsa" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.subplots(1, 1, figsize=[6, 5])\n", "ax = plt.subplot(1, 1, 1)\n", "pmsa.sel(window=1).plot()\n", "plt.title(\"PMSA\")" ] } ], "metadata": { "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.11" } }, "nbformat": 4, "nbformat_minor": 4 }