Geostrophic adjustment and inertia-gravity waves.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geostrophic adjustment and inertia-gravity waves\n",
    "### MAQ - 32806, Chiel van Heerwaarden and Imme Benedict, 2016\n",
    "In this tutorial you will use the previously derived shallow water equations to study the concept of geostrophic adjustment (again pioneered by Carl-Gustav Rossby). Geostrophic adjustment is the process that restores geostrophic and thermal wind balance after a disturbance, such as a big cloud system or a jet instability. In this tutorial, you will study how geostrophic adjustment leads to the formation of a low pressure system out of a disturbance, and you will learn about the role of inertia-gravity waves therein.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading the required packages\n",
    "First, we load the required packages and set the physical constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML\n",
    "figsize_s = 6,4\n",
    "figsize_l = 9,4\n",
    "pl.rcParams.update({'legend.fontsize' : 'medium'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Omega = 7.29e-5 # Angular velocity of the earth.\n",
    "g = 9.81        # Gravity constant.\n",
    "a = 6.371e6     # Earth radius."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### 1. The shallow water equations for small perturbations\n",
    "In this assignment, you are going to work with the linearized shallow water equations. This means that we assume that there is no background velocity, and that the velocities and height perturbations are small. Furthermore, to simplify the dynamics, we assume the Coriolis parameter $f$ to be a constant $f_0$. In addition to the derivation in class, we assume now that we have a flat bottom, with $z_b = 0$. Therefore, we can write $h$ directly into the conservation equations for horizontal momentum.\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\nonumber\n",
    "\\dfrac{\\partial u^\\prime}{\\partial t} & = & f_0 v^\\prime - g\\dfrac{\\partial h^\\prime}{\\partial x}\n",
    "\\\\\n",
    "\\nonumber\n",
    "\\dfrac{\\partial v^\\prime}{\\partial t} & = & -f_0 u^\\prime - g\\dfrac{\\partial h^\\prime}{\\partial y}\n",
    "\\\\\n",
    "\\nonumber\n",
    "\\dfrac{\\partial h^\\prime}{\\partial t} & = &\n",
    "- \\overline{h} \\left( \\dfrac{\\partial u^\\prime}{\\partial x} +\\dfrac{\\partial v^\\prime}{\\partial y} \\right)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "1. Explain the assumptions behind the equations above. Discuss in this explanation why can we neglect the advection terms if the perturbations are small, and why is $h$ replaced by $\\overline{h}$ in the right-hand side of the third equation.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. The shallow water equations for our experiment\n",
    "In our experiment, you are going to study the behavior of a two-dimensional system. We study an initially non-moving water layer in which we introduce a dip in the height that only varies with $x$. The governing equations are\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\nonumber\n",
    "\\dfrac{\\partial u^\\prime}{\\partial t} & = & f_0 v^\\prime - g \\dfrac{\\partial h^\\prime}{\\partial x}\n",
    "\\\\\n",
    "\\nonumber\n",
    "\\dfrac{\\partial v^\\prime}{\\partial t} & = & -f_0 u^\\prime\n",
    "\\\\\n",
    "\\nonumber\n",
    "\\dfrac{\\partial h^\\prime}{\\partial t} & = &\n",
    "- \\overline{h} \\dfrac{\\partial u^\\prime}{\\partial x}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "1. Explain why the governing equations can be written as above if the initial perturbation has no gradients in $y$.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Setting up the experiment\n",
    "Below, you can specify the setup of the experiment. The `f_multiplier` is initially set to `1` and can be adjusted later to make the earth spin faster or slower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f_multiplier = 1.\n",
    "lat = 45.\n",
    "f0 = f_multiplier * 2.*Omega*np.sin(np.pi/180.*lat)\n",
    "L = a * (2.*np.pi) * np.cos(np.pi/180.*lat)\n",
    "L_dist = 1.e6\n",
    "h_bar = 200.\n",
    "h_prime = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the conditions above, we can create our domain (ranging from `-L/2` to `L/2`). The initial condition of $u^{\\prime}$, $v^{\\prime}$ and $h^{\\prime}$ are stored in variables `u0`, `v0` and `h0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Here we initialize the grid and the horizontal velocities.\n",
    "nx = 1501\n",
    "x = np.linspace(-L/2., L/2., nx)\n",
    "dx = x[1]-x[0]\n",
    "u0 = np.zeros(nx)\n",
    "v0 = np.zeros(nx)\n",
    "\n",
    "# The perturbation in h has the shape of a Gaussian (normal distribution).\n",
    "x_sigma = L_dist/6.\n",
    "h0 = np.exp(-x**2/(2.*x_sigma**2))\n",
    "h0 = - h_prime*h0/np.max(h0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This plot shows the shape of the initial perturbation\n",
    "pl.figure(figsize=figsize_s)\n",
    "pl.plot(x, h0, 'b-')\n",
    "pl.xlabel('x (m)')\n",
    "pl.ylabel('h\\' (m)')\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Study the plot of the $h^{\\prime}$ as a function of $x$. Is the perturbation really small? Compare its magnitude with that of $\\overline{h}$.\n",
    "2. Calculate the initial potential vorticity $q$ as a variable `q0` and plot it as a function of $x$.\n",
    "3. Calculate the Rossby radius of deformation $\\lambda_R \\equiv \\dfrac{\\sqrt{g \\overline{h}}}{f_0}$ for this experiment.\n",
    "4. Replot $q$ and $h^{\\prime}$ with the horizontal axis normalized by the Rossby radius $L_r$ and normalize the axis between -8 and 8 $x/\\lambda_R$. How large is the perturbation compared to the Rossby radius?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Running the experiment\n",
    "Now, you will run the model for `nt` timesteps with time increments of `dt` seconds. The output is saved every `ni` steps. The final result of $u^{\\prime}$, $v^{\\prime}$ and $h^{\\prime}$ are stored in variables `u`, `v` and `h`. Below you find an animation that shows the time evolution that can aid you in interpreting what is happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nt = 50000\n",
    "ni = 500\n",
    "dt = 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# MODEL INTEGRATOR, THIS MIGHT TAKE A WHILE!\n",
    "u = np.copy(u0)\n",
    "v = np.copy(v0)\n",
    "h = np.copy(h0)\n",
    "\n",
    "h_out = np.zeros((nt//ni+1, nx))\n",
    "h_out[0,:] = h0\n",
    "\n",
    "for n in range(1,nt+1):\n",
    "    u += dt * (   f0*v - g*np.gradient(h,dx) )\n",
    "    v += dt * ( - f0*u )\n",
    "    h += dt * ( - h_bar*np.gradient(u,dx) )\n",
    "    \n",
    "    if (n%ni == 0):\n",
    "        h_out[n//ni,:] = h[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# ANIMATION CREATION PHASE 1. PLEASE IGNORE THIS FIGURE, WE DO NOT KNOW HOW TO SUPPRESS.\n",
    "# Set up the figure, the axis, and the plot element we want to animate\n",
    "fig, ax = pl.subplots(figsize=figsize_l)\n",
    "xmin, xmax = min(x), max(x)\n",
    "hmin, hmax = np.min(h_out), np.max(h_out)\n",
    "hmin -= 0.05*(hmax-hmin)\n",
    "hmax += 0.05*(hmax-hmin)\n",
    "ax.set_xlim(xmin, xmax)\n",
    "ax.set_ylim(hmin, hmax)\n",
    "ax.set_xlabel('x (m)')\n",
    "ax.set_ylabel('h (m)')\n",
    "ax.plot(x, h_out[0,:], 'b:')\n",
    "line, = ax.plot([], [], lw=2)\n",
    "fig.tight_layout()\n",
    "pl.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# ANIMATION CREATION PHASE 2\n",
    "# Initialization function: plot the background of each frame\n",
    "def init():\n",
    "    line.set_data([], [])\n",
    "    return (line,)\n",
    "\n",
    "def animate(i):\n",
    "    line.set_data(x, h_out[i,:])\n",
    "    return (line,)\n",
    "\n",
    "# Call the animator. Keyword blit=True means only re-draw the parts that have changed.\n",
    "anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                               frames=h_out.shape[0], interval=40, blit=True)\n",
    "\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### 4. Analysis of the experiment\n",
    "Below you find the end result of the experiment. All three variables $h^{\\prime}$, $u^{\\prime}$ and $v^{\\prime}$ are plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=figsize_l)\n",
    "pl.subplot(121)\n",
    "pl.plot(x, h , 'b-')\n",
    "pl.plot(x, h0, 'b:')\n",
    "pl.xlabel('x (m)')\n",
    "pl.ylabel('h\\' (m)')\n",
    "pl.subplot(122)\n",
    "pl.plot(x, u, 'b-', label='u')\n",
    "pl.plot(x, v, 'r-', label='v')\n",
    "pl.xlabel('x (m)')\n",
    "pl.ylabel('u, v (m/s)')\n",
    "pl.legend(loc=0, frameon=False);\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. What do you see in the figure above? Is this what you expected? Why does the surface not become flat?\n",
    "2. Replot the figure above, but normalize $x$ again with $\\lambda_R$ and limit the range between -8 and 8. Within how many Rossby radii is the initial perturbation contained?\n",
    "3. Calculate the potential vorticity at the end of the simulation and compare it to the initial value. How has it changed?\n",
    "4. Rerun the experiment for an earth that spins five times faster (set `f_multiplier=5.`) and for one that spins five times slower (set `f_multiplier=0.2`). What do you observe?\n",
    "5. What is the influence of the rotation rate on the phase velocity and group velocity of the waves? Can you link your observations to the dispersion relation of inertia-gravity waves?\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Geostrophic adjustment\n",
    "With the help of the potential vorticity, it is possible for this experiment to predict the final steady state out of the initial conditions. The formulation of the steady state is this specific experiment is\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\nonumber\n",
    "0 & = & f_0 v^\\prime - g \\dfrac{\\partial h^\\prime}{\\partial x}\n",
    "\\\\\n",
    "\\nonumber\n",
    "0 & = & -f_0 u^\\prime\n",
    "\\\\\n",
    "\\nonumber\n",
    "0 & = &\n",
    "- \\overline{h} \\dfrac{\\partial u^\\prime}{\\partial x}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "We know that throughout the geostrophic adjustment the potential vorticity has not changed, thus\n",
    "\n",
    "$$\n",
    "q = q_0\n",
    "$$\n",
    "___\n",
    "1. Why is the process that you have modelled called geostrophic adjustment? TIPS: Think about how to define geostrophic wind speed in the shallow water equations. Don't forget that $\\dfrac{\\partial h^{\\prime}}{\\partial y} = 0$ in this experiment.\n",
    "2. Derive the equation that relates the final height $h^{\\prime}$ to the initial potential vorticity $q_0$. This is a differential equation, you do not have to solve it.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below you find a solver for the final state as a function of the initial potential vorticity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# MAKE SURE THAT YOU HAVE CALCULATED THE INITIAL POTENTIAL VORTICITY\n",
    "# AS A VARIABLE q0\n",
    "q0 = f0 / (h_bar + h0)\n",
    "\n",
    "def calc_steady_state():\n",
    "    c1 = g/(f0*dx**2)\n",
    "\n",
    "    # Matrix solver\n",
    "    M = np.zeros((nx, nx))\n",
    "\n",
    "    for i in range(1,nx-1):\n",
    "        M[i,i-1] = c1\n",
    "        M[i,i  ] = -(2.*c1 + q0[i])\n",
    "        M[i,i+1] = c1\n",
    "\n",
    "    M[ 0, 0] = 1.\n",
    "    M[-1,-1] = 1.\n",
    "\n",
    "    h_mat = -1.*np.ones(nx)*f0 + q0*h_bar\n",
    "    h_mat[ 0] = 0\n",
    "    h_mat[-1] = 0\n",
    "\n",
    "    h = np.linalg.inv(M).dot(h_mat)\n",
    "    v = g/f0*np.gradient(h, dx)\n",
    "    \n",
    "    return h, v\n",
    "\n",
    "h_inv, v_inv = calc_steady_state()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=figsize_l)\n",
    "pl.subplot(121)\n",
    "# ADD HERE THE SOLUTION FOR H (h_inv)\n",
    "pl.plot(x, h, 'k--', label='model')\n",
    "pl.xlabel('x (m)')\n",
    "pl.ylabel('h\\' (m)')\n",
    "pl.legend(loc=0, frameon=False)\n",
    "pl.subplot(122)\n",
    "# ADD HERE THE SOLUTION FOR V (v_inv)\n",
    "pl.plot(x, v, 'k--', label='model')\n",
    "pl.xlabel('x (m)')\n",
    "pl.ylabel('v (m/s)')\n",
    "pl.legend(loc=0, frameon=False);\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Add the solution acquired from the matrix inversion to the plot above, and normalize $x$ with the Rossby radius. Is the result the same? If not, explain the differences.\n",
    "2. Explain the pattern in the final solution for $v$. What are the similarities between this pattern and a zonal cross section through the center of a low pressure system in the lower troposphere? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "End of tutorial.\n",
    "___"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}