permafrost_alaska.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "522c0574",
   "metadata": {},
   "source": [
    "# Data Component Use Case for Permafrost Thaw and Hillslope Diffusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f6caae7",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Permafrost is defined as any material (rock or soil) that remains below 0°C for two or more consecutive years. This Jupyter notebook demonstrates how to use several [CSDMS Data Components](https://csdms.colorado.edu/wiki/DataComponents) to download topography, snow, and temperature data to calculate the permafrost active layer thickness and simulate the hillslope evolution process for a study area in Alaska. \n",
    "\n",
    "In this notebook, it includes the following sections:\n",
    "- [Initial Setup](#setup)\n",
    "  \n",
    "  Install API key files and create the input/output folders.\n",
    "  <br>\n",
    "- [Step 1: Download Datasets](#step1) \n",
    "\n",
    "  Download the topography, temperature and snow datasets.\n",
    "  <br>\n",
    "- [Step 2: Calculate Activate Layer Thickness ](#step2)\n",
    "\n",
    "  Use temperature and snow datasets to calculate the activate layer thickness of the study area.  \n",
    "  <br>\n",
    "- [Step 3: Simulate Hillslope Diffusion](#step3)\n",
    "\n",
    "  Use the activate layer thickness and topography data as the input to simulate the hillslope diffusion process.\n",
    "  <br>\n",
    "\n",
    "\n",
    "**Suggested Citation**: Gan, T., Tucker, G. E., Overeem, I., Pierce, E. (2024). Data Component Use Case for Permafrost Thaw and Hillslope Diffusion, HydroShare, https://www.hydroshare.org/resource/3e11df71c1724df18fe8fed9e0afdfa1/\n",
    "\n",
    "**Run this notebook**: Please follow the instructions [here](https://github.com/gantian127/permafrost_usecase#data-component-use-case-for-permafrost-thaw-and-hillslope-diffusion) to run this notebook on the local PC or the online platform."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a9f425",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "## Initial Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9f8d9bd",
   "metadata": {},
   "source": [
    "### Install API key files\n",
    "For the ERA5 and Topography data components, there is a need to create API key files to download the datasets. The install_api_key( ) function will ask for your [CDS API Key](https://cds.climate.copernicus.eu/api-how-to) and [Open Topography API Key](https://opentopography.org/blog/introducing-api-keys-access-opentopography-global-datasets) to create API key files. Please make sure you have already obtained those API Keys before you run this helper function. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ddd865",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import install_api_key\n",
    "\n",
    "install_api_key()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f2bfe5",
   "metadata": {},
   "source": [
    "### Create folders\n",
    "We will first import all the python packages and then create three folders for this notebook:\n",
    "- **configuration file folder**: this folder includes several configuration files which will be used by the data components. In this example, we have prepared these configuration files ('dem_config.yaml' and 'era5_config.yaml') and put them in this folder. \n",
    "- **cache folder**: this folder stores the downloaded data files.\n",
    "- **results folder**: this folder stores the final results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1f07eb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages\n",
    "import os\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray\n",
    "import cftime\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "from matplotlib.colors import LightSource\n",
    "import imageio.v2 as imageio\n",
    "from IPython.display import Video\n",
    "from tqdm import trange\n",
    "\n",
    "from pymt.models import Topography, Era5, Ku\n",
    "from landlab import RasterModelGrid, imshow_grid\n",
    "from landlab.components import DepthDependentDiffuser\n",
    "\n",
    "warnings.simplefilter(action=\"ignore\", category=FutureWarning)\n",
    "plt.rcParams.update({\"font.size\": 12})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82f450b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create folders\n",
    "study_area = \"alaska\"\n",
    "\n",
    "config_dir = os.path.join(os.getcwd(), \"config_files_{}\".format(study_area))\n",
    "results_dir = os.path.join(os.getcwd(), \"results_{}\".format(study_area))\n",
    "cache_dir = os.path.join(os.getcwd(), \"cache_{}\".format(study_area))\n",
    "\n",
    "\n",
    "for folder in [config_dir, results_dir, cache_dir]:\n",
    "    if not os.path.isdir(folder):\n",
    "        os.mkdir(folder)\n",
    "        print(folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16eca92d",
   "metadata": {},
   "source": [
    "<a id='step1'></a>\n",
    "## Step 1 Download Datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f60dfec",
   "metadata": {},
   "source": [
    "### Background \n",
    "\n",
    "Permafrost covers nearly 85% of Alaska. The map below shows the permafrost distribution in this area (map source from [Jorgenson et al.](https://permafrost.gi.alaska.edu/sites/default/files/AlaskaPermafrostMap_Front_Dec2008_Jorgenson_etal_2008.pdf)). A warming climate brought higher temperatures which may cause some permafrost to thaw. This can lead to geologic hazards such as landslides, ground subsidence, erosion and other severe surface distortions.\n",
    "\n",
    "In this use case, we will run a model to calculate the active layer thickness for a study area in Alaska. Active layer is the top layer of the permafrost where the annual maximum temperature reaches 0°C and the temperature shift over diurnal and seasonal cycles. Then we will use the active layer thickness as the input to simulate the hillslope diffusion process for the study area. And we need to prepare the following datasets for this use case.\n",
    "- OpenTopography DEM \n",
    "- ERA5 temperature \n",
    "- ERA5 snow \n",
    "\n",
    "\n",
    "<img src=\"https://github.com/gantian127/permafrost_usecase/blob/master/permafrost.png?raw=true\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "573e5e1d",
   "metadata": {},
   "source": [
    "### OpenTopography DEM \n",
    "\n",
    "We will use the [Topography data component](https://csdms.colorado.edu/wiki/Model:Topography_Data_Component) to download the Digital Elevation Model (DEM) data with 30m resolution. The 'dem_config.yaml' file includes the parameter settings of this data component. The following cells demonstrate how to use the configuration file to initialize a data component and how to use the variable and grid related methods of this data component to get the metadata as well as the data values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a264880b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize Topography data component\n",
    "dem = Topography()\n",
    "dem.initialize(os.path.join(config_dir, \"dem_config.yaml\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "239ce76e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM variable info\n",
    "var_name = dem.output_var_names[0]\n",
    "var_unit = dem.var_units(var_name)\n",
    "var_location = dem.var_location(var_name)\n",
    "var_type = dem.var_type(var_name)\n",
    "var_grid = dem.var_grid(var_name)\n",
    "var_itemsize = dem.var_itemsize(var_name)\n",
    "var_nbytes = dem.var_nbytes(var_name)\n",
    "print(\n",
    "    \"variable_name: {} \\nvar_unit: {} \\nvar_location: {} \\nvar_type: {} \\nvar_grid: {} \\nvar_itemsize: {}\"\n",
    "    \"\\nvar_nbytes: {} \\n\".format(\n",
    "        var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985d5bea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM grid info\n",
    "dem_grid_ndim = dem.grid_ndim(var_grid)\n",
    "dem_grid_type = dem.grid_type(var_grid)\n",
    "dem_grid_shape = dem.grid_shape(var_grid)\n",
    "dem_grid_spacing = dem.grid_spacing(var_grid)\n",
    "dem_grid_origin = dem.grid_origin(var_grid)\n",
    "\n",
    "print(\n",
    "    \"grid_ndim: {} \\ngrid_type: {} \\ngrid_shape: {} \\ngrid_spacing: {} \\ngrid_origin: {}\".format(\n",
    "        dem_grid_ndim, dem_grid_type, dem_grid_shape, dem_grid_spacing, dem_grid_origin\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a9816f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM variable data\n",
    "dem_data = dem.get_value(var_name)\n",
    "dem_data_2D = dem_data.reshape(dem_grid_shape)\n",
    "\n",
    "# get X, Y extent for plot\n",
    "min_y, min_x = dem_grid_origin\n",
    "max_y = min_y + dem_grid_spacing[0] * (dem_grid_shape[0] - 1)\n",
    "max_x = min_x + dem_grid_spacing[1] * (dem_grid_shape[1] - 1)\n",
    "dy = dem_grid_spacing[0] / 2\n",
    "dx = dem_grid_spacing[1] / 2\n",
    "dem_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]\n",
    "\n",
    "# plot DEM data as hillshade plot\n",
    "fig, ax = plt.subplots(figsize=(8, 9))\n",
    "ls = LightSource(azdeg=90, altdeg=30)\n",
    "shaded_dem = ls.hillshade(np.array(dem_data_2D))\n",
    "\n",
    "ax.imshow(shaded_dem, cmap=\"Greys\", extent=dem_extent)\n",
    "ax.ticklabel_format(useOffset=False)\n",
    "ax.set_xlabel(\"latitude [degree_north]\")\n",
    "ax.set_ylabel(\"longitude [degree_east]\")\n",
    "ax.title.set_text(\"Study area of Eightmile Lake\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d34b5f63",
   "metadata": {},
   "source": [
    "### ERA5 Snow and Temperature\n",
    "We will use the [ERA5 data component](https://csdms.colorado.edu/wiki/Model:ERA5_Data_Component) to download the monthly mean snow water equivalent, snow density and 2m air temperature datasets for the study area with 0.25 degrees (27-28km) resolution. We will download two datasets for 1980-1989 and 2010-2019 separately.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c184492",
   "metadata": {},
   "source": [
    "'era5_1980_config.yaml' and 'era5_2010_config.yaml' files include the parameter settings for two ERA5 data components. The following cells demonstrate how to use the configuration file to initialize an ERA5 data component and how to use the variable, grid and time related methods to get the metadata as well as the data values. \n",
    "\n",
    "You'll notice that although the ERA5 and Topography data components download the datasets from different sources, they are using the same methods to get information from the datasets. Please note that sometimes the request for ERA5 data may be queued which may take a while (>10min) to get the data downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c730ff9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize ERA5 data components\n",
    "\n",
    "# 1980-1989 data\n",
    "era5 = Era5()\n",
    "era5.initialize(os.path.join(config_dir, \"era5_1980_config.yaml\"))\n",
    "\n",
    "# 2010-2019 data\n",
    "era5_2 = Era5()\n",
    "era5_2.initialize(os.path.join(config_dir, \"era5_2010_config.yaml\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cece7801",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 variable info\n",
    "for var_name in era5.output_var_names:\n",
    "    var_unit = era5.var_units(var_name)\n",
    "    var_location = era5.var_location(var_name)\n",
    "    var_type = era5.var_type(var_name)\n",
    "    var_grid = era5.var_grid(var_name)\n",
    "    var_itemsize = era5.var_itemsize(var_name)\n",
    "    var_nbytes = era5.var_nbytes(var_name)\n",
    "    print(\n",
    "        \"variable_name: {} \\nvar_unit: {} \\nvar_location: {} \\nvar_type: {} \\nvar_grid: {} \\nvar_itemsize: {}\"\n",
    "        \"\\nvar_nbytes: {} \\n\".format(\n",
    "            var_name,\n",
    "            var_unit,\n",
    "            var_location,\n",
    "            var_type,\n",
    "            var_grid,\n",
    "            var_itemsize,\n",
    "            var_nbytes,\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bc1e2e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 grid info\n",
    "era5_grid_ndim = era5.grid_ndim(var_grid)\n",
    "era5_grid_type = era5.grid_type(var_grid)\n",
    "era5_grid_shape = era5.grid_shape(var_grid)\n",
    "era5_grid_spacing = era5.grid_spacing(var_grid)\n",
    "era5_grid_origin = era5.grid_origin(var_grid)\n",
    "\n",
    "print(\n",
    "    \"grid_ndim: {} \\ngrid_type: {} \\ngrid_shape: {} \\ngrid_spacing: {} \\ngrid_origin: {}\".format(\n",
    "        era5_grid_ndim,\n",
    "        era5_grid_type,\n",
    "        era5_grid_shape,\n",
    "        era5_grid_spacing,\n",
    "        era5_grid_origin,\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "894f48bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 time info\n",
    "era5_start_time = era5.start_time\n",
    "era5_end_time = era5.end_time\n",
    "era5_time_unit = era5.time_units\n",
    "\n",
    "print(\n",
    "    \"start_time:{} \\nend_time:{} \\ntime_unit:{}\".format(\n",
    "        era5_start_time, era5_end_time, era5_time_unit\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f5583aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 variables data and plot (at the first time step)\n",
    "fig = plt.figure(figsize=(14, 12))\n",
    "nrows, ncols = 2, 2\n",
    "i = 1\n",
    "\n",
    "for var_name in era5.output_var_names:\n",
    "    ax = fig.add_subplot(nrows, ncols, i)\n",
    "    var_unit = era5.var_units(var_name)\n",
    "\n",
    "    # get variable data\n",
    "    era5_data = era5.get_value(var_name)\n",
    "    era5_data_2D = era5_data.reshape(era5_grid_shape)\n",
    "\n",
    "    # get X, Y extent for plot\n",
    "    min_y, min_x = era5_grid_origin\n",
    "    max_y = min_y + era5_grid_spacing[0] * (era5_grid_shape[0] - 1)\n",
    "    max_x = min_x + era5_grid_spacing[1] * (era5_grid_shape[1] - 1)\n",
    "    dy = era5_grid_spacing[0] / 2\n",
    "    dx = era5_grid_spacing[1] / 2\n",
    "    era5_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]\n",
    "\n",
    "    # plot data\n",
    "    im = ax.imshow(era5_data_2D, extent=era5_extent, cmap=\"Blues\")\n",
    "    ax.title.set_text(\"{} ({})\".format(var_name, var_unit))\n",
    "    ax.set_xlabel(\"longitude [degree_east]\")\n",
    "    ax.set_ylabel(\"latitude [degree_north]\")\n",
    "    cbar = plt.colorbar(im, ax=ax)\n",
    "\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6222e476",
   "metadata": {},
   "source": [
    "<a id='step2'></a>\n",
    "## Step 2 Calculate Activate Layer Thickness"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7337961",
   "metadata": {},
   "source": [
    "In this section, we will use the ERA5 datasets to prepare inputs and run Ku model (Kudryavtsev et al.,1974) to calculate the active layer thickness. The Ku model provides a steady-state solution under the assumption of sinusoidal air temperature forcing. It considers snow, vegetation, and soil layers as thermal damping to variation of air temperature. The layer of soil is considered to be a homogeneous column with different thermal properties in the frozen and thawed states. One of its main outputs is the annual active layer thickness at the top of permafrost.\n",
    "\n",
    "Given the fact that the air temperature has an increasing trend over the years, we want to explore whether this change has an impact on the active layer thickness for the study area. So we will prepare two sets of inputs for 1980-1989 and 2010-2019 separately and compare the model results.\n",
    "\n",
    "This section includes the following tasks:\n",
    "- Get ERA5 data for study area\n",
    "- Prepare temperature and snow inputs \n",
    "- Run Ku model\n",
    "- Results analysis\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec85eddb",
   "metadata": {},
   "source": [
    "### Get ERA5 data for study area\n",
    "\n",
    "From the Topography and ERA5 data plots, it can be found that the study area falls in the upper left grid of the ERA5 data. So we will get the time series of air temperature, snow density and snow depth from the two ERA5 data components for the study area. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34622d6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create dataframe to store time series data\n",
    "era5_df = pd.DataFrame(columns=[\"temp\", \"swe\", \"dens\", \"time\"])\n",
    "time_steps = 12 * 10  # 10 years of monthly data\n",
    "\n",
    "for data_comp in [era5, era5_2]:\n",
    "\n",
    "    for i in range(0, time_steps):\n",
    "        # get values\n",
    "        temp = data_comp.get_value(\"2 metre temperature\")\n",
    "        swe = data_comp.get_value(\"Snow depth\")\n",
    "        dens = data_comp.get_value(\"Snow density\")\n",
    "        time = cftime.num2pydate(data_comp.time, data_comp.time_units)\n",
    "\n",
    "        # add new row to dataframe\n",
    "        era5_df.loc[len(era5_df)] = [temp[0], swe[0], dens[0], time]\n",
    "\n",
    "        # update to next time step\n",
    "        data_comp.update()\n",
    "\n",
    "era5_df = era5_df.set_index(\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aadbafec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot time series data\n",
    "for start, end in [(1980, 1989), (2010, 2019)]:\n",
    "    ax = era5_df[(start <= era5_df.index.year) & (era5_df.index.year <= end)].plot(\n",
    "        y=[\"temp\", \"swe\", \"dens\"],\n",
    "        subplots=True,\n",
    "        figsize=(8, 7),\n",
    "        xlabel=\"\",\n",
    "        legend=None,\n",
    "        title=\"ERA5 monthly mean data for {}-{}\".format(start, end),\n",
    "    )\n",
    "    ax[0].set_title(\"air temperature (K)\")\n",
    "    ax[1].set_title(\"snow water equivalent (m)\")\n",
    "    ax[2].set_title(\"snow density (kg/m3)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24e1e63e",
   "metadata": {},
   "source": [
    "### Prepare temperature and snow inputs\n",
    "\n",
    "Now, we will use the time series data for the study area to prepare inputs for Ku model, which includes annual mean temperature, temperature amplitude and snow thickness. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcd43bd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a dataframe to store input data\n",
    "input_df = pd.DataFrame(columns=[\"temp_mean\", \"temp_amp\", \"snow_h\"])\n",
    "\n",
    "# calculate annual mean temperature and temperature amplitude\n",
    "input_df[\"temp_mean\"] = era5_df[\"temp\"].groupby(pd.Grouper(freq=\"Y\")).mean() - 273.15\n",
    "input_df[\"temp_amp\"] = (\n",
    "    era5_df[\"temp\"].groupby(pd.Grouper(freq=\"Y\")).max()\n",
    "    - era5_df[\"temp\"].groupby(pd.Grouper(freq=\"Y\")).mean()\n",
    ")\n",
    "\n",
    "# calculate snow thickness (SWE * water_dens/snow_dens)\n",
    "water_dens = 1000  # kg/m3\n",
    "snow_df = era5_df[\"swe\"] * water_dens / era5_df[\"dens\"]\n",
    "snow_df = snow_df[\n",
    "    (snow_df.index.month <= 5) | (snow_df.index.month >= 9)\n",
    "]  # remove summer month\n",
    "input_df[\"snow_h\"] = snow_df.groupby(pd.Grouper(freq=\"Y\")).mean()\n",
    "\n",
    "input_df = input_df.round(3).dropna()\n",
    "\n",
    "# plot the input datasets\n",
    "for start, end in [(1980, 1989), (2010, 2019)]:\n",
    "    ax = input_df[(start <= input_df.index.year) & (input_df.index.year <= end)].plot(\n",
    "        y=[\"temp_mean\", \"temp_amp\", \"snow_h\"],\n",
    "        subplots=True,\n",
    "        figsize=(10, 8),\n",
    "        xlabel=\"\",\n",
    "        legend=None,\n",
    "        title=\"Temperature and snow input for {}-{}\".format(start, end),\n",
    "    )\n",
    "    ax[0].set_title(\"annual mean temperature (°C)\")\n",
    "    ax[1].set_title(\"temperature amplitude (°C)\")\n",
    "    ax[2].set_title(\"snow thickness (m)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97b3cc5b",
   "metadata": {},
   "source": [
    "Let's look at these plots above. What do you find from them? It can be seen that the annual mean temperature tends to increase while the temperature amplitude and snow thickness decreased during 2010-2019. How will this change impact the active layer thickness? Will it become thicker because of temperature increase? Let's run Ku model to find it out."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85fd17bc",
   "metadata": {},
   "source": [
    "### Run Ku model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8eeaf08",
   "metadata": {},
   "source": [
    "Ku model has been implemented as one of the [pymt Model Components](https://pymt.readthedocs.io/en/latest/models.html#model-components). And we will import Ku model component from pymt for model run.\n",
    "\n",
    "Since the model or data is wrapped with the [Basic Model Interface](https://bmi.readthedocs.io/en/latest/) to become a component under the [pymt](https://pymt.readthedocs.io/en/latest/index.html) modeling framework, the way to control or query a model or a data component is the same. For example, the methods to initialize, update and get value for Ku model component are the same as the ERA5 data component. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a530779",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a dataframe to store results\n",
    "active_layer = pd.DataFrame(columns=[\"active_h\"], index=input_df.index)\n",
    "\n",
    "# run Ku model\n",
    "for start, end in [(1980, 1989), (2010, 2019)]:\n",
    "    print(\"Simulation for {}-{}\".format(start, end))\n",
    "\n",
    "    # get input data\n",
    "    input_data = input_df[(input_df.index.year >= start) & (input_df.index.year <= end)]\n",
    "\n",
    "    # setup model\n",
    "    ku = Ku()\n",
    "    args = ku.setup(start_year=start, end_year=end, lat=63.88, lon=-149.25)\n",
    "    ku.initialize(*args)\n",
    "\n",
    "    # run model\n",
    "    for index, row in input_data.iterrows():\n",
    "        ku.set_value(\"atmosphere_bottom_air__temperature\", row[\"temp_mean\"])\n",
    "        ku.set_value(\"atmosphere_bottom_air__temperature_amplitude\", row[\"temp_amp\"])\n",
    "        ku.set_value(\"snowpack__depth\", row[\"snow_h\"])\n",
    "        ku.update()\n",
    "\n",
    "        # store result\n",
    "        active_layer.loc[index] = ku.get_value(\"soil__active_layer_thickness\")[0]\n",
    "\n",
    "print(\"Simulation is done!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae6dc230",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot Ku model result\n",
    "for start, end in [(1980, 1989), (2010, 2019)]:\n",
    "    ax = active_layer[\n",
    "        (start <= active_layer.index.year) & (active_layer.index.year <= end)\n",
    "    ].plot.bar(y=[\"active_h\"], subplots=True, figsize=(7, 5), xlabel=\"\", legend=None)\n",
    "    ax[0].set_ylabel(\"Thickness (m)\")\n",
    "    ax[0].set_xticklabels(np.arange(start, end + 1, step=1), rotation=0)\n",
    "    ax[0].set_title(\"Active Layer Thickness for {}-{}\".format(start, end))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24cb9b8d",
   "metadata": {},
   "source": [
    "### Results analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21f0605a",
   "metadata": {},
   "source": [
    "In the active layer thickness plot, the change of temperature and snow conditions during 2010-2019 doesn't impact much on the active layer thickness to make it increase as we expected. \n",
    "\n",
    "Why is that? To find out the reason, we will first calculate the 10-year average of the input and output data for 1980-1989 and 2010-2019 and do some model experiments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12067c1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get 10 year average of input and output for 1980-1989 and 2010-2019\n",
    "ave_df = pd.DataFrame(\n",
    "    columns=[\"temp_mean\", \"temp_amp\", \"snow_h\", \"active_h\"], index=[1980, 2010]\n",
    ")\n",
    "\n",
    "for start, end in [(1980, 1989), (2010, 2019)]:\n",
    "    ave_input = input_df[\n",
    "        (input_df.index.year >= start) & (input_df.index.year <= end)\n",
    "    ].mean()\n",
    "    ave_output = active_layer[\n",
    "        (start <= active_layer.index.year) & (active_layer.index.year <= end)\n",
    "    ].mean()\n",
    "    ave_df.loc[start] = [\n",
    "        ave_input[\"temp_mean\"],\n",
    "        ave_input[\"temp_amp\"],\n",
    "        ave_input[\"snow_h\"],\n",
    "        ave_output[\"active_h\"],\n",
    "    ]\n",
    "\n",
    "ave_df.astype(\"float\").round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32d2207b",
   "metadata": {},
   "source": [
    "The 10-year average data proves what we found from the temperature and snow inputs plot and the active layer thickness plot. The increased annual mean temperature and decreased temperature amplitude as well as snow thickness didn't lead to much change of the active layer thickness (about 4mm difference). \n",
    "\n",
    "We will run several model experiments to find out the impact of those model inputs on the active layer thickness. We will run the model using 10-year average input for 1980-1989 and 2010-2019. The model run for 1980-1989 will be taken as the \"base\" experiment for comparison. Based on the model inputs of 1980-1989, we will change its annual mean temperature, temperature amplitude and snow thickness value separately by replacing them with the 10-year average input of 2010-2019."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1cb1305",
   "metadata": {},
   "outputs": [],
   "source": [
    "# experiment inputs\n",
    "experiment = {\n",
    "    \"1980-1989\": [-3.849, 16.99, 0.355],\n",
    "    \"2010-2019\": [-2.563, 16.234, 0.247],\n",
    "    \"change of annual mean temperature\": [-2.563, 16.99, 0.355],\n",
    "    \"change of temperature amplitude\": [-3.849, 16.234, 0.355],\n",
    "    \"change of snow thickness\": [-3.849, 16.99, 0.247],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0274eb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run experiments\n",
    "for key, value in experiment.items():\n",
    "    ku = Ku()\n",
    "    args = ku.setup(\n",
    "        T_air=value[0], A_air=value[1], h_snow=value[2], lat=63.88, lon=-149.25\n",
    "    )\n",
    "    ku.initialize(*args)\n",
    "    ku.update()\n",
    "\n",
    "    print(\n",
    "        \"Result for {}: {}m\".format(\n",
    "            key, ku.get_value(\"soil__active_layer_thickness\")[0].round(3)\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d64bb4c",
   "metadata": {},
   "source": [
    "The experiment results showed that it can lead to an increase of the active layer thickness by only increasing the annual temperature. But especially if the snow thickness decreases, its insulating capacity in mid and late winter will decrease, and as such the active layer thickness will also decrease. Therefore, the respective change of warming temperature versus a decreasing snow thickness can act in opposing direction and thereby minimize changes for the active layer thickness. This phenomenon was also observed with field datasets and studied by several researchers at other study sites (Garnello et al., 2021; Zhang, 2005).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96e9588c",
   "metadata": {},
   "source": [
    "<a id='step3'></a>\n",
    "## Step 3  Simulate Hillslope Diffusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26df085c",
   "metadata": {},
   "source": [
    "In this section, we will use the topography DEM and the active layer thickness as the inputs for a hillslope diffusion model. The [Landlab](https://landlab.readthedocs.io/en/master/) component [DepthDependentDiffuser]() will be used for simulation. This component implements a depth and slope dependent linear diffusion rule in the style of Johnstone and Hilley (2015). \n",
    "\n",
    "This section includes the following tasks:\n",
    "- Create model grid and data fields\n",
    "- Initialize and run model component \n",
    "- Results visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393e88e1",
   "metadata": {},
   "source": [
    "### Create model grid and data fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "675f6f8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a model grid\n",
    "model_grid = RasterModelGrid(dem_grid_shape, xy_spacing=30)\n",
    "\n",
    "# add soil depth field and set it as active layer thickness\n",
    "soil_depth = model_grid.add_zeros(\"soil__depth\", at=\"node\", clobber=True)\n",
    "\n",
    "active_h = 0.992  # result for 2010-2019 with 10-year average input\n",
    "model_grid.at_node[\"soil__depth\"].fill(active_h)\n",
    "\n",
    "# add soil production rate field and set it as 0\n",
    "production = model_grid.add_zeros(\"soil_production__rate\", at=\"node\", clobber=True)\n",
    "\n",
    "# add elevation field using topography DEM data\n",
    "elevation = model_grid.add_field(\n",
    "    \"node\", \"topographic__elevation\", dem_data.astype(\"float\"), units=\"m\", clobber=True\n",
    ")\n",
    "\n",
    "# plot elevation data field\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "ax = plt.gca()\n",
    "ax.invert_yaxis()\n",
    "imshow_grid(\n",
    "    model_grid,\n",
    "    \"topographic__elevation\",\n",
    "    plot_name=\"Elevation of Eight Mile Lake\",\n",
    "    var_name=\"elevation (m)\",\n",
    ")\n",
    "ax.set_xlabel(\"east-west distance (m)\")\n",
    "ax.set_ylabel(\"north-south distance (m)\")\n",
    "fig.savefig(os.path.join(results_dir, \"elev_0.png\".format(time)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc0ddd35",
   "metadata": {},
   "source": [
    "### Initialize and run model component "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2040ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set parameters\n",
    "linear_diffusivity = 1e-2  # m^2/yr\n",
    "transport_decay_depth = 1  # m\n",
    "dt = 1  # time step (year)\n",
    "nt = 1000  # total steps (years)\n",
    "\n",
    "# get original dem data\n",
    "ref_data = np.copy(elevation)\n",
    "\n",
    "# initialize the landlab component\n",
    "ddiff = DepthDependentDiffuser(model_grid, linear_diffusivity, transport_decay_depth)\n",
    "\n",
    "for time in trange(1, nt + 1):\n",
    "    # run model\n",
    "    ddiff.run_one_step(dt)\n",
    "\n",
    "    # save result\n",
    "    if time % 50 == 0:\n",
    "        im_data = model_grid.at_node[\"topographic__elevation\"] - ref_data\n",
    "\n",
    "        fig = plt.figure(figsize=(8, 6))\n",
    "        ax = plt.gca()\n",
    "        ax.invert_yaxis()\n",
    "        imshow_grid(\n",
    "            model_grid,\n",
    "            im_data,\n",
    "            plot_name=\"Simulated elevation change after {} years\".format(time),\n",
    "            var_name=\"elevation (m)\",\n",
    "            vmin=-0.1,\n",
    "            vmax=0.1,\n",
    "            cmap=\"winter_r\",\n",
    "        )\n",
    "        ax.set_xlabel(\"east-west distance (m)\")\n",
    "        ax.set_ylabel(\"north-south distance (m)\")\n",
    "        plt.close(fig)\n",
    "        fig.savefig(os.path.join(results_dir, \"elev_{}.png\".format(time)))\n",
    "\n",
    "print(\"Simulation is done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cffe17bc",
   "metadata": {},
   "source": [
    "### Results visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09af4607",
   "metadata": {},
   "source": [
    "Let's run the cells bellow to visualize the results. The plot shows the simulated elevation after 1000 years of hillslope diffusion process. The video shows the elevation change during the simulation period. You will find the vertical stripe pattern from the video which is mainly caused by the similar pattern from the terrain elevation data and can be improved with other datasets in a higher resolution (e.g., Lidar DEM dataset). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f1f5903",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# plot result\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "ax = plt.gca()\n",
    "ax.invert_yaxis()\n",
    "imshow_grid(\n",
    "    model_grid,\n",
    "    \"topographic__elevation\",\n",
    "    plot_name=\"Simulated Elevation after 1000 years\",\n",
    "    var_name=\"elevation (m)\",\n",
    ")\n",
    "ax.set_xlabel(\"east-west distance (m)\")\n",
    "ax.set_ylabel(\"north-south distance (m)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d7a2372",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a short video\n",
    "img_files = [\n",
    "    os.path.join(results_dir, file)\n",
    "    for file in os.listdir(results_dir)\n",
    "    if \".png\" in file\n",
    "]\n",
    "img_files.sort(key=lambda x: os.path.getmtime(x))\n",
    "\n",
    "with imageio.get_writer(\n",
    "    os.path.join(results_dir, \"hillslope.mp4\"), mode=\"I\", fps=1, macro_block_size=None\n",
    ") as writer:\n",
    "    for f in img_files:\n",
    "        im = imageio.imread(os.path.join(results_dir, f))\n",
    "        writer.append_data(im)\n",
    "\n",
    "writer.close()\n",
    "\n",
    "# Display the video\n",
    "Video(\"./results_alaska/hillslope.mp4\", embed=True, width=500, height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2af04b5b",
   "metadata": {},
   "source": [
    "## References\n",
    "- Anisimov, O. A., Shiklomanov, N. I., & Nelson, F. E. (1997). Global warming and active-layer thickness: results from transient general circulation models. Global and Planetary Change, 15(3-4), 61-77. https://doi.org/10.1016/S0921-8181(97)00009-X\n",
    "- Garnello, A., Marchenko, S., Nicolsky, D., Romanovsky, V., Ledman, J., Celis, G., Schädel, C., Luo, Y., & Schuur, E. A. G. (2021). Projecting Permafrost Thaw of Sub-Arctic Tundra With a Thermodynamic Model Calibrated to Site Measurements. Journal of Geophysical Research: Biogeosciences, 126(6), e2020JG006218. https://doi.org/10.1029/2020JG006218\n",
    "- Johnstone, S., Hilley, G. (2015). Lithologic control on the form of soil-mantled hillslopes, Geology 43(1), 83-86. https://doi.org/10.1130/G36052.1\n",
    "- Kudryavtsev, V. A., Garagulya, L. S., Kondrat'yeva, K. A., Melamed, V. G. (1974). Fundamentals of frost forecasting in geological engineering investigations Nauka, Moscow, p. 431. https://apps.dtic.mil/sti/pdfs/ADA039677.pdf\n",
    "- Zhang, T. (2005). Influence of the seasonal snow cover on the ground thermal regime: An overview. Reviews of Geophysics, 43(4). https://doi.org/10.1029/2004RG000157\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30fdf5b8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}