Exercise-5-problems-1-3.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "5001ba7f00883c8ed72315cfe72d91d4",
          "grade": false,
          "grade_id": "cell-f91fdb5001d60d49",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "3UmA5uZXuBpr"
      },
      "source": [
        "# Exercise 5: Viscous flow of ice (20 points)\n",
        "\n",
        "This week's exercise explores the flow of ice through a valley as a channel and down a slope.\n",
        "\n",
        "### Tips for completing this exercise\n",
        "\n",
        "- Use **exactly** the same variable names as in the instructions because your answers will be automatically graded, and the tests that grade your answers rely on following the same formatting or variable naming as in the instructions.\n",
        "- **Please do not**:\n",
        "\n",
        "    - **Change the file names**. Do all of your editing in the provided `Exercise-5-problems-1-3.ipynb` file (this file).\n",
        "    - **Copy/paste cells in this notebook**. We use an automated grading system that will fail if there are copies of code cells.\n",
        "    - **Change the existing cell types**. You can add cells, but changing the cell types for existing cells (from code to markdown, for example) will also cause the automated grader to fail."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "2ae8159beb12ad38c3763c4ee9af1067",
          "grade": false,
          "grade_id": "cell-005b6534a2f43b1b",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "VvuPLwCsuBpz"
      },
      "source": [
        "## Problem 1 - Ice flow in an open channel (12 points)\n",
        "\n",
        "The goal of this problem is to calculate horizontal velocity profiles across a glacier for Newtonian and non-Newtonian fluid flow resulting from a pressure gradient.\n",
        "\n",
        "**Your score on this problem will be based on following criteria**:\n",
        "\n",
        "- Creating a Python function for calculating the pressure gradient\n",
        "- Creating a Python function for the channel flow velocity\n",
        "- Calculating and plotting ice velocities in a channel\n",
        "- Including comments that explain what most lines in the code do\n",
        "- Uploading your notebook and `introqg_functions.py` script file to your GitHub repository for this week's exercise"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "2eba8ecde61bd95837ca20467ec97597",
          "grade": false,
          "grade_id": "cell-9f80085b4268c165",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "LAHegUQsuBp1"
      },
      "source": [
        "### Part 0: Copying and testing your script file from Exercise 4 (0 points)\n",
        "\n",
        "The first task in this problem is to copy your `introqg_functions.py` script file from Exercise 4 to the directory containing this notebook and then run the cell below to ensure it has been copied and is functioning as expected. Note: We will only check some of the functions in this file using the tests below, not all of them.\n",
        "\n",
        "- Copy your `introqg_functions.py` script file from Exercise 4 to the directory containing this notebook\n",
        "- Run the tests below"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "4090d0f97662e91ce3dc31675a6b372d",
          "grade": false,
          "grade_id": "cell-6480d1978df821d2",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "Vav0BGdouBp2"
      },
      "outputs": [],
      "source": [
        "# The test below should work\n",
        "\n",
        "from nose.tools import assert_equal\n",
        "from introqg_functions import steady_state_temp, transient_temp\n",
        "\n",
        "# Test the steady-state temperature calculation function\n",
        "test_temp1 = steady_state_temp(\n",
        "    gradient=10.0,\n",
        "    diffusivity=32.0,\n",
        "    velocity=0.5,\n",
        "    depths=10.0)\n",
        "\n",
        "test_temp2 = steady_state_temp(\n",
        "    gradient=20.0,\n",
        "    diffusivity=32.0,\n",
        "    velocity=0.75,\n",
        "    depths=15.0)\n",
        "\n",
        "test_temp3 = transient_temp(\n",
        "    initial_gradient=20.0,\n",
        "    diffusivity=32.0,\n",
        "    velocity=0.5,\n",
        "    depths=10.0,\n",
        "    time=5.0)\n",
        "\n",
        "test_temp4 = transient_temp(\n",
        "    initial_gradient=15.0,\n",
        "    diffusivity=32.0,\n",
        "    velocity=1.5,\n",
        "    depths=15.0,\n",
        "    time=15.0)\n",
        "\n",
        "# Print calculated closure temperature\n",
        "print(f\"Calculated steady-state temperature 1: {test_temp1:7.3f} °C\")\n",
        "print(f\"Calculated steady-state temperature 2: {test_temp2:7.3f} °C\")\n",
        "print(f\"Calculated transient temperature 1: {test_temp3:7.3f} °C\")\n",
        "print(f\"Calculated transient temperature 2: {test_temp4:7.3f} °C\")\n",
        "\n",
        "# Check that the closure temperature values are correct\n",
        "assert_equal(round(test_temp1, 3), 92.579)\n",
        "assert_equal(round(test_temp2, 3), 252.938)\n",
        "assert_equal(round(test_temp3, 3), 232.717)\n",
        "assert_equal(round(test_temp4, 3), 465.720)\n",
        "\n",
        "# Print message if it is safe to continue\n",
        "print(\"\\nAll tests pass! You are ready to proceed with this exercise.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "1cad4175e96418ac69cdaa4bed87f345",
          "grade": false,
          "grade_id": "cell-0713deadb7cdb9f9",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "s3J0pYM9uBp6"
      },
      "source": [
        "### Part 1: Background and getting started (0 points)\n",
        "\n",
        "Before you start you should read through the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-1) for this problem. The other parts of this problem will refer to items described in this theory text.\n",
        "\n",
        "- Read through the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-1) for this problem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "0bcd4da77c14ad5b3b7a9b3ebb38ced3",
          "grade": false,
          "grade_id": "cell-936f5cea4af18b4c",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "J-gPOFEAuBp7"
      },
      "source": [
        "### Part 2: A pressure gradient function? (4 points)\n",
        "\n",
        "As mentioned in the description of the problem, we will be calculating velocities of ice at the surface as we cross a glacial valley. We can calculate the ice velocity using Equation 10 from the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-1), which is also shown below\n",
        "\n",
        "\\begin{equation}\n",
        "  \\large\n",
        "  u(y) = \\frac{a}{n + 1} \\left( \\frac{p_{1} - p_{0}}{L} \\right)^{n} \\left[ \\left( \\frac{h}{2} \\right)^{n + 1} - y^{n + 1} \\right]\n",
        "\\end{equation}\n",
        "\n",
        "where $a$ is a function of the viscosity and temperature (we will ignore temperature for now), $n$ is the power-law exponent, $\\frac{p_{1} - p_{0}}{L}$ is the pressure gradient in the channel, $h$ is the channel width, and $y$ is array for distance from the center of the channel. The only problem here is that the value for the pressure gradient $\\left(\\frac{p_{1} - p_{0}}{L}\\right)^{n}$ is not known. Thus, our first step is to create a function to estimate its value from some of the observed velocities for the Saskatchewan glacier.\n",
        "\n",
        "- Create a function in your `introqg_functions.py` script file called `channel_pressure` that calculates the value of $\\left(\\frac{p_{1} - p_{0}}{L}\\right)^{n}$ for a given value of the power-law exponent `n`.\n",
        "    - You can solve for $\\left(\\frac{p_{1} - p_{0}}{L}\\right)^{n}$ by rearranging Equation 10 (see above) to solve for $\\left(\\frac{p_{1} - p_{0}}{L}\\right)^{n}$ and using the condition that the maximum flow velocity $u = 2.41 \\times 10^{-6}$ m/s occurs at $y = 0$ (see Table 1 in Part 4).\n",
        "    - Assume $a = 5 \\times 10^{-24}$ Pa<sup>-3</sup> s<sup>-1</sup> and $h = 1400$ m.\n",
        "    - Your function parameters should include\n",
        "        - The power-law exponent `exponent`\n",
        "        - The channel width `width`\n",
        "        - The viscosity/temperature constant `viscosity_constant`\n",
        "        - The channel velocity at y = 0 `velocity_max`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "1f57d2fe0f734342d9160a2d05a7d1f6",
          "grade": false,
          "grade_id": "problem_2_ag_p2_pressure",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "sMz_h6tQuBp9"
      },
      "outputs": [],
      "source": [
        "# Import your channel_pressure function below\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "39cf77e023efd061bdc56f4dc05956f0",
          "grade": true,
          "grade_id": "problem_2_ag_p2_pressure_test1",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "dRgNqR4yuBp-"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "from nose.tools import assert_equal\n",
        "\n",
        "# Function test calculations\n",
        "pressure1 = channel_pressure(\n",
        "    exponent=3,\n",
        "    width=1400.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=2.41e-6\n",
        ")\n",
        "\n",
        "# Print results\n",
        "print(f\"The pressure gradient for exponent=3 is {pressure1:.3f}. Expected value: 8029987.505.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(pressure1, 3), 8029987.505)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "30c773311736e152392bf3e5feb045e2",
          "grade": true,
          "grade_id": "problem_2_ag_p2_pressure_test2",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "A1eeNIlUuBp_"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "\n",
        "# Function test calculations\n",
        "pressure2 = channel_pressure(\n",
        "    exponent=6,\n",
        "    width=1400.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=2.41e-6\n",
        ")\n",
        "pressure3 = channel_pressure(\n",
        "    exponent=1.5,\n",
        "    width=1400.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=2.41e-6\n",
        ")\n",
        "\n",
        "# Print results\n",
        "print(f\"The pressure gradient for exponent=6 is {pressure2:.3f}. Expected value: 0.041.\")\n",
        "print(f\"The pressure gradient for exponent=1.5 is {pressure3:.3f}. Expected value: 92948406117.575.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(pressure2, 3), 0.041)\n",
        "assert_equal(round(pressure3, 3), 92948406117.575)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "8d4e6536ce1517a5231b11c9c8b67193",
          "grade": false,
          "grade_id": "cell-22aef2bdb4215115",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "-2_F9J5guBqA"
      },
      "source": [
        "### Part 3: An ice flow velocity function (4 points)\n",
        "\n",
        "Your next task is to create a function to calculate the velocity of ice flowing in a channel (Equation 10).\n",
        "\n",
        "- Create a function in your `introqg_functions.py` script file called `channel_velocity` that calculates the velocity of ice flowing in a channel using Equation 10 from the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-1), which was shown in the previous part.\n",
        "    - Your function should use the `channel_pressure` function created in the previous part.\n",
        "    - As in the previous part you should assume $a = 5 \\times 10^{-24}$ Pa<sup>-3</sup> s<sup>-1</sup>, $h = 1400$ m, and that the flow is symmetric about $y = 0$.\n",
        "    - Your function parameters should include\n",
        "        - The power-law exponent `exponent`\n",
        "        - The channel width array `y`\n",
        "        - The channel width `width`\n",
        "        - The viscosity/temperature constant `viscosity_constant`\n",
        "        - The channel velocity at y = 0 `velocity_max`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "302042509427a6dc0428c43becc64404",
          "grade": false,
          "grade_id": "problem_2_ag_p3_velocity",
          "locked": false,
          "schema_version": 3,
          "solution": true
        },
        "id": "A-FCgQkquBqB"
      },
      "outputs": [],
      "source": [
        "# Import your channel_velocity function below\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "991ba6aa920deed3926d545e505cae7e",
          "grade": true,
          "grade_id": "problem_2_ag_p3_velocity_test1",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false
        },
        "id": "5rVMRwjIuBqC"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "import numpy as np\n",
        "from nose.tools import assert_equal\n",
        "\n",
        "# Fake data\n",
        "y1 = np.array([-500.0, -250.0, 0.0, 200.0, 400.0, 600.0])\n",
        "velocity1 = channel_velocity(exponent=2, y=y1)\n",
        "\n",
        "# Print results\n",
        "print(f\"The channel velocity at y1[1] for exponent=2 is {velocity1[1]:.5e} m/s. Expected value: 2.30022e-06 m/s.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(velocity1[-1], 8), 8.9e-07)\n",
        "assert_equal(round(velocity1[3], 8), 2.35e-06)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "7119b0f9cea0984401f68853674c33dc",
          "grade": true,
          "grade_id": "problem_2_ag_p3_velocity_test2",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "aA6SYsQjuBqD"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "\n",
        "# Fake data\n",
        "y2 = np.array([-300, -150.0, 0.0, 100.0, 200.0, 300.0])\n",
        "velocity2 = channel_velocity(exponent=4, y=y2)\n",
        "\n",
        "# Print results\n",
        "print(f\"The channel velocity at y2[1] for exponent=4 is {velocity2[1]:.5e} m/s. Expected value: 2.40891e-06 m/s.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(velocity2[-1], 8), 2.38e-06)\n",
        "assert_equal(round(velocity2[3], 8), 2.41e-06)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "9de73fbdcd760f7532f397aa30c46857",
          "grade": false,
          "grade_id": "cell-a0bb7a59ed8c5aef",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "WgFVjGgAuBqD"
      },
      "source": [
        "### Part 4: Calculating ice flow velocities (2 points)\n",
        "\n",
        "Your next task is to calculate the flow velocity for ice flowing in a channel using your `channel_velocity` function.\n",
        "\n",
        "- Create a new pandas DataFrame called `channel` to store channel velocities for 1400-m-wide ice channel\n",
        "    - Use `y` as the index for your `channel` DataFrame\n",
        "    - **Hint**: You should define $width = 1400$ m and assume that the flow is symmetric about $y = 0$ and define `y` with a length of 201\n",
        "- Calculate the velocity of ice in a channel as a non-Newtonian fluid with power-law exponents of $exponent$ = 2, 3, 4, and 5.\n",
        "    - As above, you can assume $a = 5 \\times 10^{-24}$ Pa<sup>-3</sup> s<sup>-1</sup>\n",
        "- Store the output from each call of the `channel_velocity` function for each value of $exponent$ in your `channel` DataFrame\n",
        "    - You should label the new column in your DataFrame using the power-law exponent. For example, you could have a column called `exponent = 2` for a power-law exponent of 2.\n",
        "    - **Note**: You should convert your channel velocities from m/s to m/yr before storing them in the `channel` DataFrame!\n",
        "    - The hints for this week's exercise may be helpful for this task"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "1539dfe6c8acd104e0ca488524ef3442",
          "grade": false,
          "grade_id": "problem_2_ag_p4_velocities",
          "locked": false,
          "schema_version": 3,
          "solution": true
        },
        "id": "T3UivjQUuBqE"
      },
      "outputs": [],
      "source": [
        "# Import pandas\n",
        "import pandas as pd\n",
        "\n",
        "# Enter your values below\n",
        "width = None\n",
        "y = None\n",
        "channel = None\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "aa5e085706940358e96ff9ddda924199",
          "grade": true,
          "grade_id": "problem_2_ag_p4_velocities_test1",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false
        },
        "id": "4-9qtnoCuBqE"
      },
      "outputs": [],
      "source": [
        "# These tests should work\n",
        "print(f\"The channel shape is {channel.shape}. Expected value: (201, 4).\")\n",
        "print(f\"The velocity at channel.iloc[51,0] is {channel.iloc[51,0]:.4f} m/yr. Expected value: 67.1062 m/yr.\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "cd32f448cb79fb179837e0e71ca2c1ef",
          "grade": false,
          "grade_id": "cell-025b16de6a269a19",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "0y3cxhN6uBqF"
      },
      "source": [
        "### Part 5: Plotting your results (2 points)\n",
        "\n",
        "Next you can create a plot of the velocity calculations you've done.\n",
        "\n",
        "- Create a plot of the velocity of the ice for all 4 power-law exponents\n",
        "    - Be sure to label your axes and add a title.\n",
        "    - Also be sure your line legend makes it clear which profile is which\n",
        "- Add a figure caption in the Markdown cell beneath the plot, describing the figure as if it were in a scientific journal article"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "7f77f81833891a58b2e9b53c520cd697",
          "grade": true,
          "grade_id": "problem_2_mg_p5_plot",
          "locked": false,
          "points": 1.5,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "souKkxw-uBqG"
      },
      "outputs": [],
      "source": [
        "# Add your plot below\n",
        "ax = None\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "c04d96fe7be4019dc14f6bfbaddb0036",
          "grade": true,
          "grade_id": "problem_2_mg_p5_plot_caption",
          "locked": false,
          "points": 0.5,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "wIx5C7ZiuBqH"
      },
      "source": [
        "YOUR ANSWER HERE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "6107ab7613a7c1ee1ea99d52ebec689e",
          "grade": false,
          "grade_id": "cell-85bc31e63965a084",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "bf0O1wH4uBqI"
      },
      "source": [
        "## Problem 2 - Comparing our model to the Saskatchewan glacier velocities (8 points)\n",
        "\n",
        "In this problem we will use the `channel_velocity` function from Problem 2 to compare predicted ice velocities to glacier flow velocities observed from the [Saskatchewan glacier near Banff in Alberta, Canada](https://goo.gl/maps/R9S48J4KYPr) (Figure 1). The Saskatchewan glacier is 1400 m wide and part of a large ice field known as the Columbia Icefield.\n",
        "\n",
        "**Your score on this problem will be based on following criteria**:\n",
        "\n",
        "- Reading in the data file\n",
        "- Comparing predictions for the ice flow velocity to observations from the data file\n",
        "- Plotting your results\n",
        "- Including comments that explain what most lines in the code do\n",
        "- Uploading your notebook and `introqg_functions.py` script file to your GitHub repository for this week's exercise\n",
        "\n",
        "![Saskatchewan glacier](https://upload.wikimedia.org/wikipedia/commons/a/a5/Saskatchewan_Glacier.jpg)<br/>\n",
        "*Figure 1. Saskatchewan glacier.*."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "c8e28389b47f2a72b6cd9a06daab0383",
          "grade": false,
          "grade_id": "cell-8fe062ca0d4f7d9a",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "O5eVPai9uBqI"
      },
      "source": [
        "### Part 1: Reading the data file (2 points)\n",
        "\n",
        "Distance (m) | Velocity (m/yr) | Velocity (m/s)\n",
        "------------ | -------------- | --------------\n",
        "-660 | 12 | 3.80E-07\n",
        "-640 | 22 | 7.00E-07\n",
        "-570 | 42 | 1.33E-06\n",
        "-460 | 63 | 2.00E-06\n",
        "-220 | 74 | 2.35E-06\n",
        "40 | 76 | 2.41E-06\n",
        "180 | 74 | 2.35E-06\n",
        "260 | 72 | 2.28E-06\n",
        "500 | 51 | 1.62E-06\n",
        "\n",
        "*Table 1. Velocity measurements across the Saskatchewan glacier.*\n",
        "\n",
        "The file [`data/sask_glacier_velo.txt`](data/sask_glacier_velo.txt) contains a series of surface velocities measured at various locations across the glacier (Table 1). Your task now is to load in this velocity data.\n",
        "\n",
        "- Read the data file [`data/sask_glacier_velo.txt`](data/sask_glacier_velo.txt) using pandas into a DataFrame called `sask`\n",
        "    - Set the column `'Distance (m)'` as the index column"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "ccd7cc79aea0397eeca93f3081af2ff6",
          "grade": false,
          "grade_id": "problem_3_ag_p2_read",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "IcGXYNixuBqJ"
      },
      "outputs": [],
      "source": [
        "# Load the data file\n",
        "sask = None\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "6fde48583cf47b6b442c832bd945e99a",
          "grade": true,
          "grade_id": "problem_3_ag_p2_read_test1",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "QEJQDZ_XuBqJ"
      },
      "outputs": [],
      "source": [
        "# This test should work\n",
        "print(f\"The first 5 rows of the data file are:\\n\\n{sask.head()}\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "134807862d56d11b2aeca89005773b51",
          "grade": false,
          "grade_id": "cell-8fe062ca0d4f7d9ab",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "QT8z_fDDuBqK"
      },
      "source": [
        "### Part 2: Comparing to observations (4 points)\n",
        "\n",
        "The next task is to calculate ice velocities at locations corresponding to the observed velocity measurements and compare the predictions to the data in the data file.\n",
        "\n",
        "- Loop over the power-law exponent values $exponent=2$ to $exponent=5$ and calculate the velocity of ice using the `channel_velocity` function at the positions where the velocity of the Saskatchewan glacier has been observed in the data file\n",
        "    - Use the default function values other than the power-law exponent and values for `y`\n",
        "    - You should label the new column in your `sask` DataFrame using the power-law exponent. For example, you could have a column called `exponent = 2` for a power-law exponent of 2.\n",
        "    - **Note**: You should convert your channel velocities from m/s to m/yr before storing them in the `sask` DataFrame!\n",
        "- For each set of predicted ice flow velocities, calculate the goodness-of-fit to the observed ice velocities and append this to the `gof_sask` list\n",
        "    - We will be using your `chi_squared` function from back in Exercise 2\n",
        "    - We are assuming an arbitrary standard deviation for the observed velocities of 2 m/yr.\n",
        "        - When you call your `chi_squared` function you can pass the variable `sask_stdev` for the standard deviations. It is defined in the cell below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "9685a8f994bfe3b358b3912a24d0ee55",
          "grade": false,
          "grade_id": "problem_3_ag_p3_gof",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "SVheDg-ouBqK"
      },
      "outputs": [],
      "source": [
        "# Import items we need\n",
        "import numpy as np\n",
        "from introqg_functions import chi_squared\n",
        "\n",
        "# Create array of standard deviations\n",
        "sask_stdev = np.ones(len(sask)) * 2.0\n",
        "\n",
        "# Create empty list for goodness-of-fit values\n",
        "gof_sask = []\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "f404cec5a14b5ca0d3af5cb751b29b6b",
          "grade": true,
          "grade_id": "problem_3_ag_p3_gof_test1",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "KPz0RGh5uBqL"
      },
      "outputs": [],
      "source": [
        "# This test should work\n",
        "print(f\"The first calculated ice velocity is {sask.iloc[0,2]:.2f} m/yr. Expected value: 12.31 m/yr.\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "0f253244b0fb9e47890ec6df4d7016ce",
          "grade": true,
          "grade_id": "problem_3_ag_p3_gof_test2",
          "locked": true,
          "points": 2,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "3KBcVd4iuBqL"
      },
      "outputs": [],
      "source": [
        "# This test should work\n",
        "print(f\"The first goodness of fit value is {gof_sask[0]:.2f}. Expected value: 4.06.\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "b7d8d793e852618d5f3151a1433e85b7",
          "grade": false,
          "grade_id": "cell-e19f29f8c7e83ed7",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "xIvOfsuIuBqL"
      },
      "source": [
        "### Part 3: Plotting your results (2 points)\n",
        "\n",
        "The last task in this problem is to plot the predicted ice flow velocities along with the observed values and show the best goodness of fit as text on the plot.\n",
        "\n",
        "- Plot the measured velocities along with the 4 predicted velocity profiles you have calculated **in meters per year (m/yr)**\n",
        "    - For the predicted velocities, use the values from the `channel` DataFrame from Problem 1\n",
        "    - For the observed velocities, use the values from the `sask` DataFrame\n",
        "        - Plot these as black circles\n",
        "    - Be sure to label your axes and include a title\n",
        "    - Also add the best goodness-of-fit value as text on the plot along with the corresponding power-law exponent value\n",
        "- Add a figure caption in the Markdown cell beneath the plot, describing the figure as if it were in a scientific journal article"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "34177727c5d73bc0b44b6729bd0902f6",
          "grade": true,
          "grade_id": "problem_3_mg_p4_plot",
          "locked": false,
          "points": 1.5,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "sAXx-wBquBqM"
      },
      "outputs": [],
      "source": [
        "# Plot your results below\n",
        "ax = None\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "a6c51fdc13677d568dd118190a721d04",
          "grade": true,
          "grade_id": "problem_3_mg_p4_plot_caption",
          "locked": false,
          "points": 0.5,
          "schema_version": 3,
          "solution": true
        },
        "id": "nnYzxsvTuBqM"
      },
      "source": [
        "YOUR ANSWER HERE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "14ff22fcce3aa082b693df47f66ba24f",
          "grade": false,
          "grade_id": "cell-728758e35292a008",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "BBMF3N_ruBqM"
      },
      "source": [
        "## Problem 3 - Non-Newtonian ice flow down an inclined plane (*Optional*, 0 points)\n",
        "\n",
        "The [Athabasca glacier](https://goo.gl/maps/HggYfoKxEUQ2) (Figure 2) is another glacier in the Columbia Icefield, which will be the focus in this problem.\n",
        "\n",
        "**Your score on this problem will be based on following criteria**:\n",
        "\n",
        "- Creating a Python function for calculating the gravitational force\n",
        "- Creating a Python function for calculating the ice velocity\n",
        "- Reading in the data file\n",
        "- Comparing predictions for the ice flow velocity to observations from the data file\n",
        "- Plotting your results\n",
        "- Including comments that explain what most lines in the code do\n",
        "- Uploading your notebook and `introqg_functions.py` script file to your GitHub repository for this week's exercise\n",
        "\n",
        "![Athabasca glacier](https://upload.wikimedia.org/wikipedia/commons/4/41/Icefields_parkway.jpg)<br/>\n",
        "*Figure 2. Athabasca glacier.*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "d5eb886e6dc9955e447c97cc6007b7f3",
          "grade": false,
          "grade_id": "cell-a7703cd4cbdc83b6",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "OthJNjUcuBqN"
      },
      "source": [
        "### Part 1: Background and theory (0 points)\n",
        "\n",
        "First read through the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-2) for this problem before proceeding.\n",
        "\n",
        "- Read through the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-2) for this problem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "49a1926af2110f534ebb832b499e0c5f",
          "grade": false,
          "grade_id": "cell-5f2e375ae6d36ab9",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "DTsdpjBzuBqN"
      },
      "source": [
        "### Part 2: A gravitational force function? (0 points)\n",
        "\n",
        "In this problem we will be calculating the flow of ice down a slope resulting from a gravitational forces on the glacier. As above, we will consider both Newtonian and non-Newtonian fluids, this time using Equation 19 from the [background and theory for this problem](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-2), which is also shown below\n",
        "\n",
        "\\begin{equation}\n",
        "  \\large\n",
        "  u(z) = u_{\\mathrm{b}} + \\frac{a}{n + 1} \\gamma_{x}^{n} \\left[ h^{n + 1} - (h - z)^{n + 1} \\right]\n",
        "\\end{equation}\n",
        "\n",
        "where $u_{\\mathrm{b}}$ is the basal sliding velocity of the glacier, $a$ is a function of the viscosity and temperature (we will ignore temperature for now), $n$ is the power-law exponent, $\\gamma_{x}$ is the downslope component of the gravitational force, $h$ is the thickness of the ice perpendicular to the slope, and $z$ is the array of distances from the base of the glacier. Similar to Problem 3, the value for the gravitational force $\\gamma_{x}$ is not known, but we can create a function to estimate its value.\n",
        "\n",
        "- Create a function in your `introqg_functions.py` script file called `gravity_force` that calculates the value of $\\gamma_{x}^{n}$ for a given value of the power-law exponent `n`.\n",
        "    - You can solve for $\\gamma_{x}^{n}$ by rearranging Equation 19 (see above) to solve for $\\gamma_{x}^{n}$ and using the condition that the maximum flow velocity $u = 9.1 \\times 10^{-7}$ m/s occurs at $z = h = 209$ m (see Table 2 in Part 4).\n",
        "    - Note that you should assume the basal sliding velocity $u_{\\mathrm{b}}$ is zero for the Newtonian case and equal to the observed sliding velocity $u_{\\mathrm{b}} = 1.30 \\times 10^{7}$ m/s for the non-Newtonian cases (see Table 2 in Part 4).\n",
        "    - Again, assume $a = 5 \\times 10^{-24}$ Pa<sup>-3</sup> s<sup>-1</sup> and $h = 209$ m.\n",
        "    - Your function parameters should include\n",
        "        - The power-law exponent `exponent`\n",
        "        - The ice thickness `thickness`\n",
        "        - The viscosity/temperature constant `viscosity_constant`\n",
        "        - The maximum ice flow velocity `velocity_max`\n",
        "        - The basal sliding velocity `velocity_basal`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "025368cd42243b555d8c973dc0d9272c",
          "grade": false,
          "grade_id": "problem_4_ag_p2_gamma_x",
          "locked": false,
          "schema_version": 3,
          "solution": true
        },
        "id": "ves62YpeuBqO"
      },
      "outputs": [],
      "source": [
        "# Import your gravity_force function below\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "a271b10c146564cd4b07c72819d20cb3",
          "grade": true,
          "grade_id": "problem_4_ag_p2_gamma_x_test1",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false
        },
        "id": "wqJXvB_wuBqO"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "from nose.tools import assert_equal\n",
        "\n",
        "# Function test calculations\n",
        "force1 = gravity_force(\n",
        "    exponent=1,\n",
        "    thickness=209.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=9.1e-7,\n",
        "    velocity_basal=1.3e-7\n",
        ")\n",
        "\n",
        "# Print results\n",
        "print(f\"The gravity force for n=1 is {force1:.2f}. Expected value: 8333142556260.16.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(force1, 2), 8333142556260.16)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "dc3273e1bd4433f12e803b7b4b1bc416",
          "grade": true,
          "grade_id": "problem_4_ag_p2_gamma_x_test2",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "xQScYZbbuBqO"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "from nose.tools import assert_equal\n",
        "\n",
        "# Function test calculations\n",
        "force2 = gravity_force(\n",
        "    exponent=5,\n",
        "    thickness=209.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=9.1e-7,\n",
        "    velocity_basal=1.3e-7\n",
        ")\n",
        "force3 = gravity_force(\n",
        "    exponent=7,\n",
        "    thickness=209.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=9.1e-7,\n",
        "    velocity_basal=1.3e-7\n",
        ")\n",
        "\n",
        "# Print results\n",
        "print(f\"The gravity force for n=5 is {force2:.2f}. Expected value: 11230.48.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(force2, 2), 11230.48)\n",
        "assert_equal(round(force3, 2), 0.34)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "2ddf9929d9b34be87ea772f7d9bd827d",
          "grade": false,
          "grade_id": "cell-83a717adf05cc85f",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "0BfMMaXRuBqP"
      },
      "source": [
        "### Part 3: An ice flow velocity function, again (0 points)\n",
        "\n",
        "Your next task is to create a function to calculate the velocity profile for ice flowing down a slope (Equation 19), shown in Part 2.\n",
        "\n",
        "- Create a function in your `introqg_functions.py` script file called `slope_velocity` that calculates the velocity of ice flowing down a slope using Equation 19 from the [background and theory](https://introqg-site.readthedocs.io/en/latest/lessons/L5/exercise-5-theory.html#problem-2), which was shown in the previous part.\n",
        "    - Your function should use the `gravity_force` function created in the previous part.\n",
        "    - As in the previous part you should assume $a = 5 \\times 10^{-24}$ Pa<sup>-3</sup> s<sup>-1</sup>, $h = 209$ m, and the maximum flow velocity $u_{\\mathrm{max}} = 9.1 \\times 10^{-7}$ m/s.\n",
        "    - Also you should again assume the basal sliding velocity $u_{\\mathrm{b}}$ is zero for the Newtonian case and equal to the observed sliding velocity $u_{\\mathrm{b}} = 1.30 \\times 10^{7}$ m/s for the non-Newtonian cases (see Table 2 in Part 4).\n",
        "    - Your function parameters should include\n",
        "        - The power-law exponent `exponent`\n",
        "        - The channel thickness array `z`\n",
        "        - The ice thickness `thickness`\n",
        "        - The viscosity/temperature constant `viscosity_constant`\n",
        "        - The maximum ice flow velocity `velocity_max`\n",
        "        - The basal sliding velocity `velocity_basal`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "520ec6b6989d65728d232348190f3b8e",
          "grade": false,
          "grade_id": "problem_4_ag_p3_slope_velocity",
          "locked": false,
          "schema_version": 3,
          "solution": true
        },
        "id": "0B-I6VUruBqP"
      },
      "outputs": [],
      "source": [
        "# Import your slope_velocity function below\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "f78b7178495d891d4a1d2445a964430c",
          "grade": true,
          "grade_id": "problem_4_ag_p3_slope_velocity_test1",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false
        },
        "id": "vSMcj64TuBqP"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "import numpy as np\n",
        "from nose.tools import assert_equal\n",
        "\n",
        "# Function test calculations\n",
        "z1 = np.array([0.0, 21.0, 74.0, 104.0, 159.0, 209.0])\n",
        "velocity1 = slope_velocity(\n",
        "    exponent=3,\n",
        "    z=z1,\n",
        "    thickness=209.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=9.1e-7,\n",
        "    velocity_basal=1.3e-7\n",
        ")\n",
        "\n",
        "# Print results\n",
        "print(f\"The velocity velocity1[4] for exponent=3 is {velocity1[4]:.8f} m/s. Expected value: 0.00000091 m/s.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(velocity1[-1], 8), 9.1e-07)\n",
        "assert_equal(round(velocity1[3], 8), 8.6e-07)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "92d6f5ca804d81b54638ef02c380fe9f",
          "grade": true,
          "grade_id": "problem_4_ag_p3_slope_velocity_test2",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "I99IgS0BuBqP"
      },
      "outputs": [],
      "source": [
        "# The visible tests below should work\n",
        "\n",
        "# Function test calculations\n",
        "velocity2 = slope_velocity(\n",
        "    exponent=1,\n",
        "    z=z1,\n",
        "    thickness=209.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=9.1e-7,\n",
        "    velocity_basal=1.3e-7\n",
        ")\n",
        "velocity3 = slope_velocity(\n",
        "    exponent=7,\n",
        "    z=z1,\n",
        "    thickness=209.0,\n",
        "    viscosity_constant=5.0e-24,\n",
        "    velocity_max=9.1e-7,\n",
        "    velocity_basal=1.3e-7\n",
        ")\n",
        "\n",
        "# Print results\n",
        "print(f\"The velocity velocity2[4] for exponent=1 is {velocity2[4]:.8f} m/s. Expected value: 0.00000086 m/s.\")\n",
        "\n",
        "# Check the other test calculations\n",
        "assert_equal(round(velocity2[-1], 8), 9.1e-07)\n",
        "assert_equal(round(velocity2[3], 8), 6.8e-07)\n",
        "assert_equal(round(velocity3[-1], 9), 9.1e-07)\n",
        "assert_equal(round(velocity3[3], 9), 9.07e-07)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "c58e2e07cc4bfa912604acf9e73d66e2",
          "grade": false,
          "grade_id": "cell-3dc006a66003d55f",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "gceIGvODuBqQ"
      },
      "source": [
        "### Part 4: Reading the data file (0 points)\n",
        "\n",
        "Depth from surface (m) | Height from base (m) | Velocity (m/yr) | Velocity (m/s)\n",
        "---------------------- | -------------------- | -------------- | --------------\n",
        "0 | 209 | 28.6 | 9.10E-07\n",
        "15 | 195 | 28.5 | 9.00E-07\n",
        "30 | 180 | 28.5 | 9.00E-07\n",
        "45 | 165 | 28.4 | 9.00E-07\n",
        "60 | 150 | 28.2 | 8.90E-07\n",
        "75 | 135 | 28.0 | 8.90E-07\n",
        "90 | 120 | 27.7 | 8.80E-07\n",
        "105 | 105 | 27.2 | 8.60E-07\n",
        "120 | 90 | 26.5 | 8.40E-07\n",
        "135 | 75 | 25.5 | 8.10E-07\n",
        "150 | 60 | 24.0 | 7.60E-07\n",
        "165 | 45 | 21.5 | 6.80E-07\n",
        "180 | 30 | 17.5 | 5.50E-07\n",
        "195 | 15 | 10 | 3.20E-07\n",
        "209 | 0 | 4 | 1.30E-07\n",
        "\n",
        "*Table 2. Velocities measured from a vertical profile through Athabasca glacier.*\n",
        "\n",
        "A vertical velocity profile for the Athabasca glacier has been measured and the measurements are in the file [`data/atha_glacier_velo.txt`](data/atha_glacier_velo.txt) (Table 2). Your task now is to load in this velocity data.\n",
        "\n",
        "- Read the data file [`data/atha_glacier_velo.txt`](data/sask_glacier_velo.txt) using pandas into a DataFrame called `athabasca`\n",
        "    - **Note**: Do NOT assign an index column, this will cause trouble when plotting the data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "05ca7ecac0c0be0045816e48c21a29ea",
          "grade": false,
          "grade_id": "problem_4_ag_p4_data",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "5NbAvvhYuBqQ"
      },
      "outputs": [],
      "source": [
        "# Import pandas\n",
        "import pandas as pd\n",
        "\n",
        "# Read the data file below\n",
        "athabasca = None\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "820f3340a2a9ae0184694c0e9c67bc3a",
          "grade": true,
          "grade_id": "problem_4_ag_p4_data_test1",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "hW1baNWruBqQ"
      },
      "outputs": [],
      "source": [
        "# This test should work\n",
        "print(f\"The first 5 rows of the data file are\\n\\n{athabasca.head()}\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "c2d1644ec58cd322137cc29629e57796",
          "grade": false,
          "grade_id": "cell-f788ae8db415e230",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "9Z-SO-7XuBqR"
      },
      "source": [
        "### Part 5: Comparing to observations, again (0 points)\n",
        "\n",
        "The next task is to calculate ice velocities at locations corresponding to the observed velocity measurements and compare the predictions to the data in the data file.\n",
        "\n",
        "- Loop over the power-law exponent values $n=1$ to $n=5$ and calculate the velocity of ice using the `slope_velocity` function at the positions where the velocity of the Athabasca glacier has been observed in the data file\n",
        "    - Use the default function values other than the power-law exponent and values for `z`\n",
        "        - **Hint**: Use the column `'Height from base (m)'` for the values for `z`\n",
        "    - You should label the new column in your `athabasca` DataFrame using the power-law exponent. For example, you could have a column called `n = 2` for a power-law exponent of 2.\n",
        "    - **Note**: You should convert your ice velocities from m/s to m/yr before storing them in the `athabasca` DataFrame!\n",
        "- For each set of predicted ice flow velocities, calculate the goodness-of-fit to the observed ice velocities and append this to the `gof_athabasca` list\n",
        "    - We will be using your `chi_squared` function from back in Exercise 2\n",
        "    - We are assuming an arbitrary standard deviation for the observed velocities of 1 m/yr.\n",
        "        - When you call your `chi_squared` function you can pass the variable `athabasca_stdev` for the standard deviations. It is defined in the cell below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "b0f4085b3dd5a3cf9bda1989e939798e",
          "grade": false,
          "grade_id": "problem_4_ag_p5_gof",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "_RX7Lai5uBqR"
      },
      "outputs": [],
      "source": [
        "# Import items we need\n",
        "import numpy as np\n",
        "from introqg_functions import chi_squared\n",
        "\n",
        "# Create array of standard deviations\n",
        "athabasca_stdev = np.ones(len(athabasca))\n",
        "\n",
        "# Create empty list for goodness-of-fit values\n",
        "gof_athabasca = []\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "f5b6fc103b5d7ea6192507825a1e3272",
          "grade": true,
          "grade_id": "problem_4_ag_p5_gof_test1",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "S5YxYt4muBqS"
      },
      "outputs": [],
      "source": [
        "# This test should work\n",
        "print(f\"The first calculated ice velocity is {athabasca.iloc[0,4]:.2f} m/yr. Expected value: 28.72 m/yr.\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "c7975e939cb09ca82deadcc2564081de",
          "grade": true,
          "grade_id": "problem_4_ag_p5_gof_test2",
          "locked": true,
          "points": 0,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "IYCg3v4uuBqS"
      },
      "outputs": [],
      "source": [
        "# This test should work\n",
        "print(f\"The last goodness of fit is {gof_athabasca[-1]:.2f}. Expected value: 1.18.\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "77704fe9b58cb03f328395aadd216355",
          "grade": false,
          "grade_id": "cell-295d4c9cf915b366",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "CNpsSUNFuBqS"
      },
      "source": [
        "### Part 6: Plotting your results (0 points)\n",
        "\n",
        "The last task in this problem is to again plot the predicted ice flow velocities along with the observed values and show the best goodness of fit as text on the plot.\n",
        "\n",
        "- Plot the measured velocities along with the 5 predicted velocity profiles you have calculated **in meters per year (m/yr)** on the x-axis as a function of their **height from the base** on the y-axis\n",
        "    - For the predicted velocities, use the values from the `athabasca` DataFrame\n",
        "    - For the observed velocities, also use the values from the `athabasca` DataFrame\n",
        "        - Plot these as black circles connected by a black line\n",
        "    - Be sure to label your axes and include a title\n",
        "    - Also add the best goodness-of-fit value as text on the plot along with the corresponding power-law exponent value\n",
        "    - **Hint**: You should use the `label` parameter for your plots to provide the text to use for the legend\n",
        "- Add a figure caption in the Markdown cell beneath the plot, describing the figure as if it were in a scientific journal article"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "1d030ead591e7d61d7d3922769b79c79",
          "grade": true,
          "grade_id": "problem_4_mg_p6_plot",
          "locked": false,
          "points": 0,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "Anbzp0W2uBqS"
      },
      "outputs": [],
      "source": [
        "# Plot your results below\n",
        "ax = None\n",
        "\n",
        "# YOUR CODE HERE\n",
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "1ff7334e72aa50c919baaa9f7eb1ef50",
          "grade": true,
          "grade_id": "problem_4_mg_p6_plot_caption",
          "locked": false,
          "points": 0,
          "schema_version": 3,
          "solution": true
        },
        "id": "tMOmd80fuBqT"
      },
      "source": [
        "YOUR ANSWER HERE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "98a48344117c2cfda18acefe700f2a97",
          "grade": false,
          "grade_id": "cell-c3fcafc16cc7d9df",
          "locked": true,
          "schema_version": 3,
          "solution": false
        },
        "id": "EXTXu5RVuBqT"
      },
      "source": [
        "## Optional reflection questions (0 points)\n",
        "\n",
        "1. Based on your models, what is the most likely power-law exponent for glacial ice?\n",
        "2. Looking at your plots, how sensitive is the velocity of ice to the power-law exponent?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "164ea26f043c6e94a1a6e44330be21b4",
          "grade": true,
          "grade_id": "summary_questions_mg",
          "locked": false,
          "points": 0,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "E3yww4JpuBqT"
      },
      "source": [
        "YOUR ANSWER HERE"
      ]
    }
  ],
  "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.8.12"
    },
    "colab": {
      "name": "Exercise-5-problems-1-3.ipynb",
      "provenance": [],
      "toc_visible": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}