docs/thorpe_scales.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thorpe Scales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "This section is available as a jupyter notebook in the `docs/` directory or at https://github.com/modscripps/mixsea/tree/main/docs.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we'll examine some of the different parameter options used for estimating the turbulent dissipation rate, $\\epsilon$, from Thorpe scales. The method was first published by Stephen Thorpe - see <cite data-cite=\"Thorpe1977\">Thorpe (1977)</cite>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Basic idea**\n",
    "\n",
    "Thorpe's method for estimating dissipation relies on the Ozmidov scale\n",
    "    \n",
    "\\begin{equation}\n",
    "L_{O}=\\alpha \\epsilon^{1/2} N^{-3/2},\n",
    "\\end{equation}\n",
    "    \n",
    "defined by dimensional analysis as the scale at which stratification begins to affect turbulence, where $\\alpha=\\mathcal{O}(1)$ is a proportionality constant and $N$ is a measure of the background stratification (we will see later the importance of its definition).\n",
    "    \n",
    "Thorpe's insight was to identify the Ozmidov scale with the scale $L_{T}$ (now named after him) of the largest vertical overturns (gravitational instability) observed in a density profile. He proposed sorting the unstable profile to produce a stable one, keeping track of individual water parcels vertical displacements to use their root mean square displacement $d(z)$ as a proxy for the overturn scale\n",
    "\\begin{equation}\n",
    "L_{T} = \\sqrt{\\overline{d^{2}}}.\n",
    "\\end{equation}\n",
    "\n",
    "Setting $L_{O} = L_{T}$ and inverting for $\\epsilon$ gives\n",
    "\\begin{equation}\n",
    "\\boxed{\\epsilon = C_{0} L_{T}^{2} N^{3},}\n",
    "\\end{equation}\n",
    "    \n",
    "where we identify $C_{0} = \\alpha^{-2}$ as an $\\mathcal{O}(1)$ free parameter, perhaps tuned based on microstructure measurements of $\\epsilon$.\n",
    "\n",
    "---\n",
    "\n",
    "***Warning:** The Thorpe scale method is an **event-based** turbulence parameterization, so one must be wary of spatial and temporal intermittency when interpreting results.*\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Thorpe scale method in practice can be separated into three steps: 1) finding the overturns, 2) determining their Thorpe scales by sorting, and 3) determining the background stratification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "First, we load modules and example data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mixsea as mx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ctd = mx.helpers.read_ctd_testfile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our example data contain some NaN values, as is common with observational datasets. Below we isolate the good data by removing NaNs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "notnan = np.isfinite(ctd[\"depth\"]) & np.isfinite(ctd[\"t\"]) & np.isfinite(ctd[\"SP\"])\n",
    "\n",
    "depth = ctd[\"depth\"][notnan]\n",
    "t = ctd[\"t\"][notnan]\n",
    "SP = ctd[\"SP\"][notnan]\n",
    "lon = ctd[\"lon\"][0]\n",
    "lat = ctd[\"lat\"][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A detailed look at the output and diagnostics\n",
    "\n",
    "Now we run the Thorpe scale function on the good test data, selecting input parameters that are reasonable for the deep ocean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dnoise = 5e-4  # Noise parameter\n",
    "alpha = 0.95  # Coefficient relating the Thorpe and Ozmidov scales.\n",
    "# Background value of epsilon applied where no overturns are detected.\n",
    "background_eps = np.nan\n",
    "# Do not use the intermediate profile method\n",
    "use_ip = False\n",
    "# Critical value of the overturn ratio\n",
    "Roc = 0.3\n",
    "\n",
    "# Calculate Thorpe scales and diagnostics.\n",
    "eps, N2, diag = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    dnoise=dnoise,\n",
    "    alpha=alpha,\n",
    "    Roc=Roc,\n",
    "    background_eps=background_eps,\n",
    "    use_ip=use_ip,\n",
    "    return_diagnostics=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard output `eps` and `N2` (buoyancy frequency squared) have been quality controlled. The type and strictness of quality control can be set using arguments to `eps_overturn` (see the function docstring for details). In addition, numerous diagnostic variables are output in the dictionary `diag`. The diagnostic dictionary contains output from pre-quality control, as well as flags which indicate the quality control state and miscellaneous other variables that might be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for key in diag:\n",
    "    print(key)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the diagnostics are as follows:\n",
    "\n",
    "* `eps`: dissipation rate estimated before applying quality control flags\n",
    "* `N2`: buoyancy frequency squared estimated before applying quality control flags\n",
    "* `Lt`: Thorpe scale\n",
    "* `thorpe_disp`: Thorpe displacement\n",
    "* `dens`: potential density (this may appear 'steppy' due to the use of different reference pressures for different pressure ranges)\n",
    "* `dens_sorted`: sorted potential density\n",
    "* `dens_ip`: density from the intermediate profile method (does not appear in the list above because the intermediate profile method was not used)\n",
    "* `t_ip`: temperature from the intermediate profile method (does not appear in the list above because intermediate profile method was not used)\n",
    "* `Ro`: overturn ratio\n",
    "\n",
    "Some of the quality control flags:\n",
    "\n",
    "* `noise_flag`: if density/CT difference from the top to bottom of the overturn is less than the noise parameter, this flag will be True\n",
    "* `N2_flag`: if buoyancy frequency is negative, this flag will be True\n",
    "* `ends_flag`: if an overturn contains the first or last point in the data, this flag will be True\n",
    "* `Ro_flag`: if the overturn ratio is less than the critical value, this flag will be True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we zoom into a segment, it is easier to understand what these diagnostics mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot only in the depth range:\n",
    "cut = (depth > 4260) & (depth < 4500)\n",
    "zc = depth[cut]\n",
    "\n",
    "# finding unsorting indices for intuitive Thorpe displacements\n",
    "# idx_unsorted = np.argsort(idx_sorted)\n",
    "\n",
    "# axs[1].plot(thorpe_displacement[idx_unsorted][inoverturn_unsorted], z[inoverturn_unsorted])\n",
    "\n",
    "fig, axs = plt.subplots(1, 4, sharey=True, figsize=(9, 6))\n",
    "axs[0].plot(diag[\"dens\"][cut] - 1000, zc, label=\"density\", lw=1.0)\n",
    "axs[0].plot(diag[\"dens_sorted\"][cut] - 1000, zc, label=\"sorted density\", lw=1.0)\n",
    "axs[1].plot(diag[\"eps\"][cut], zc, label=\"rejected epsilon\", color=\"r\")\n",
    "axs[1].plot(eps[cut], zc, label=\"accepted epsilon\", color=\"g\")\n",
    "axs[2].plot(diag[\"N2\"][cut], zc, label=\"rejected N2\", color=\"r\")\n",
    "axs[2].plot(N2[cut], zc, label=\"accepted N2\", color=\"g\")\n",
    "axs[3].plot(diag[\"thorpe_disp\"][cut], zc, label=\"Thorpe displacement\")\n",
    "axs[3].plot(diag[\"Lt\"][cut], zc, label=\"±Thorpe scale\")\n",
    "axs[3].plot(-diag[\"Lt\"][cut], zc, \"C1--\")\n",
    "axs[0].invert_yaxis()\n",
    "axs[0].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[1].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[2].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[3].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[0].set_ylabel(\"Depth [m]\")\n",
    "axs[0].set_xlabel(\"Potential density [kg m$^{-3}$]\")\n",
    "axs[1].set_xlabel(r\"$\\epsilon$ [W kg$^{-1}$]\")\n",
    "axs[2].set_xlabel(r\"$N^2$ [rad$^2$ s$^{-2}$]\")\n",
    "axs[3].set_xlabel(\"Thorpe displacement [m]\")\n",
    "\n",
    "fig.align_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from the left, we have potential density and sorted potential density. Overturning patches should correspond to regions where the potential density is different from the sorted potential density. Next is the dissipation rate, $\\epsilon$. Plotted in green are values that have passed quality control and in red those that were rejected by some criteria. The plot looks a bit odd because there are jumps in $\\epsilon$. This is caused by small overturns adjacent to large overturns. Ultimately, a lot of small overturns will be rejected by quality control. Third from the left is buoyancy frequency, which can be calculated in a number of different ways. We will explore this in more detail later. Last are Thorpe displacements and Thorpe length scale. Large negative values atop large positive values represent large overturns; whether this leads to large dissipation depends on the stratification.\n",
    "\n",
    "The state (True/False) of some of the quality control flags is plotted below for the same segment of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(\n",
    "    1, 5, sharey=True, figsize=(9, 6), gridspec_kw={\"width_ratios\": [2, 1, 1, 1, 2]}\n",
    ")\n",
    "axs[0].plot(diag[\"thorpe_disp\"][cut], zc, label=\"Thorpe displacement\")\n",
    "axs[0].plot(diag[\"Lt\"][cut], zc, label=\"Thorpe scale\")\n",
    "axs[0].plot(-diag[\"Lt\"][cut], zc, \"C1--\")\n",
    "axs[1].plot(diag[\"noise_flag\"][cut], zc, color=\"r\")\n",
    "axs[2].plot(diag[\"N2_flag\"][cut], zc, color=\"r\")\n",
    "axs[3].plot(diag[\"ends_flag\"][cut], zc, color=\"r\")\n",
    "axs[4].plot(diag[\"Ro\"][cut], zc, label=\"Ro value\")\n",
    "axs[4].plot(diag[\"Ro_flag\"][cut], zc, color=\"r\", label=\"Ro flag\")\n",
    "axs[4].axvline(Roc, color=\"k\", ls=\":\", label=\"Ro critical\")\n",
    "axs[0].invert_yaxis()\n",
    "axs[0].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[4].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[0].set_ylabel(\"Depth [m]\")\n",
    "axs[0].set_xlabel(\"Thorpe displacement [m]\")\n",
    "axs[1].set_xlabel(\"Noise flag\")\n",
    "axs[2].set_xlabel(\"N2 flag\")\n",
    "axs[3].set_xlabel(\"Ends flag\")\n",
    "axs[4].set_xlabel(\"Ro flag and value\")\n",
    "\n",
    "fig.align_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The noise flag is raised (True) where the density difference (top to bottom) of the overturning region is smaller than the noise parameter. The N2 flag has not been raised anywhere, indicating there are no regions of negative buoyancy frequency. The ends flag is raised for the last overturn since it contains the end point of the data. However, it is just a warning flag and not used to quality control the output. The Ro flag is raised where the overturn ratio is less than the critical value <cite data-cite=\"Gargett2008\">Gargett and Garner (2008)</cite>. Note that for this demonstration we have set the critical overturn ratio a bit higher than would be normal.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The intermediate profile method\n",
    "\n",
    "The intermediate profile method attempts to reduce the influence of random noise when estimating overturns. For technical details, see <cite data-cite=\"Ferron1998\">Ferron et al. (1998)</cite> and <cite data-cite=\"Gargett2008\">Gargett and Garner (2008)</cite>.\n",
    "\n",
    "What does the method do? Here we construct three intermediate profiles from a segment of temperature data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = 1e-3\n",
    "hinge = 0\n",
    "\n",
    "# Restrict ourselves to:\n",
    "cut = (depth < 4300) & (depth > 4250)\n",
    "\n",
    "tup = mx.overturn.intermediate_profile(t[cut], acc, hinge, \"up\")\n",
    "tdo = mx.overturn.intermediate_profile(t[cut], acc, hinge, \"down\")\n",
    "tav = mx.overturn.intermediate_profile(t[cut], acc, hinge, \"ave\")\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.invert_yaxis()\n",
    "ax.plot(t[cut], depth[cut], label=\"data\")\n",
    "ax.plot(tup, depth[cut], label=\"up\")\n",
    "ax.plot(tdo, depth[cut], label=\"down\")\n",
    "ax.plot(tav, depth[cut], label=\"ave\")\n",
    "ax.legend()\n",
    "ax.set_xlabel(r\"Temperature [$^\\circ$C]\")\n",
    "ax.set_ylabel(\"Depth [m]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method creates a new profile that moves in increments of the argument `acc`, which is a measure of the accuracy of the observations. A region that looks unstable in the data (such as the cold water over light in the region around 4270 m) might be smoothed out in the intermediate profile because the fluctuations are less than the accuracy. The method can be applied top down (`down`), bottom up (`up`) or the average of up and down (`ave`). \n",
    "\n",
    "Now we investigate what the method does to estimates of $\\epsilon$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dnoise = 5e-4\n",
    "alpha = 0.95\n",
    "# Use the intermediate profile method\n",
    "use_ip = True\n",
    "\n",
    "eps, N2, diag = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    dnoise=dnoise,\n",
    "    alpha=alpha,\n",
    "    use_ip=use_ip,\n",
    "    return_diagnostics=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the same depth range as the case that did not use the method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cut = (depth > 4260) & (depth < 4500)\n",
    "zc = depth[cut]\n",
    "\n",
    "fig, axs = plt.subplots(1, 4, sharey=True, figsize=(9, 6))\n",
    "axs[0].plot(diag[\"dens\"][cut] - 1000, zc, label=\"density\")\n",
    "axs[0].plot(diag[\"dens_sorted\"][cut] - 1000, zc, label=\"sorted density\")\n",
    "axs[0].plot(diag[\"dens_ip\"][cut] - 1000, zc, label=\"intermediate density\")\n",
    "axs[1].plot(diag[\"eps\"][cut], zc, label=\"rejected epsilon\", color=\"r\")\n",
    "axs[1].plot(eps[cut], zc, label=\"accepted epsilon\", color=\"g\")\n",
    "axs[2].plot(diag[\"N2\"][cut], zc, label=\"rejected N2\", color=\"r\")\n",
    "axs[2].plot(N2[cut], zc, label=\"accepted N2\", color=\"g\")\n",
    "axs[3].plot(diag[\"thorpe_disp\"][cut], zc, label=\"Thorpe displacement\")\n",
    "axs[3].plot(diag[\"Lt\"][cut], zc, label=\"±Thorpe scale\")\n",
    "axs[3].plot(-diag[\"Lt\"][cut], zc, \"C1--\")\n",
    "axs[0].invert_yaxis()\n",
    "axs[0].legend(loc=\"lower left\", bbox_to_anchor=(-0.1, 1))\n",
    "axs[1].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[2].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[3].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[0].set_ylabel(\"Depth [m]\")\n",
    "axs[0].set_xlabel(\"Potential density [kg m$^{-3}$]\")\n",
    "axs[1].set_xlabel(r\"$\\epsilon$ [W kg$^{-1}$]\")\n",
    "axs[2].set_xlabel(r\"$N^2$ [rad$^2$ s$^{-2}$]\")\n",
    "axs[3].set_xlabel(\"Thorpe displacement [m]\")\n",
    "\n",
    "fig.align_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're using the intermediate profile method, the results have changed dramatically. Overturns are diagnosed from the intermediate profile of density, which is different to the actual density. The overturning patches may be of different sizes to those estimated without the intermediate profile. Sorting the intermediate profile results in strange looking Thorpe displacements, because the sorting algorithm struggles when many data points have an identical value. There is an argument to be made the method should not be applied to data that has already been bin averaged (as ours has) and the method may work better on raw CTD data that is more noisy.\n",
    "\n",
    "Because the intermediate profile method fundamentally changes many aspects of the Thorpe scale method, including the size of detected overturns and the Thorpe displacements, it is not straightforward or necessarily meaningful to raise a flag for rejected overturns. Consequently, most flags are disabled when using the method and so are not shown."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare methods of estimating $N^2$\n",
    "\n",
    "There are a few different options for calculating buoyancy frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_ip = False\n",
    "\n",
    "# Choose buoyancy frequency method\n",
    "N2_method = \"teosp1\"\n",
    "eps_teosp1, N2_teosp1 = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    ")\n",
    "\n",
    "N2_method = \"teos\"\n",
    "eps_teos, N2_teos = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    dnoise,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    ")\n",
    "\n",
    "N2_method = \"bulk\"\n",
    "eps_bulk, N2_bulk = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    ")\n",
    "\n",
    "N2_method = \"endpt\"\n",
    "eps_endpt, N2_endpt = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A brief explanation of the methods:\n",
    "\n",
    "* `teos`: This method estimates $N^2$ from the thermodynamic equation of state 2010 (TEOS), using the tempurature and salinity at the first and last point in the overturn.\n",
    "* `teosp1`: This method is the same as above except it uses the temperature and salinity from the points immediately above and below the overturn (p1 meaning plus 1).\n",
    "* `bulk`: This method comes from <cite data-cite=\"Smyth2001\">Smyth et al. (2001)</cite> and has the advantage of being insensitive to errors in determining the patch size. Buoyancy frequency is proportional to the root mean square density anomaly of an overturn, divided by the Thorpe scale.\n",
    "* `endpt`: This simple method estimates buoyancy frequency from the potential density gradient across each overturn.\n",
    "\n",
    "More methods may be added in future versions of `mixsea`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot only in the depth range:\n",
    "cut = (depth > 4225) & (depth < 4350)\n",
    "zc = depth[cut]\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, sharey=True, figsize=(9, 6))\n",
    "axs[0].plot(diag[\"dens\"][cut] - 1000, zc, label=\"density\")\n",
    "axs[0].plot(diag[\"dens_sorted\"][cut] - 1000, zc, label=\"sorted density\")\n",
    "axs[1].plot(eps_teosp1[cut], zc, label=\"teosp1\")\n",
    "axs[1].plot(eps_teos[cut], zc, label=\"teos\")\n",
    "axs[1].plot(eps_bulk[cut], zc, label=\"bulk\")\n",
    "axs[1].plot(eps_endpt[cut], zc, label=\"endpt\")\n",
    "axs[2].plot(N2_teosp1[cut], zc, label=\"teosp1\")\n",
    "axs[2].plot(N2_teos[cut], zc, label=\"teos\")\n",
    "axs[2].plot(N2_bulk[cut], zc, label=\"bulk\")\n",
    "axs[2].plot(N2_endpt[cut], zc, label=\"endpt\")\n",
    "axs[0].invert_yaxis()\n",
    "axs[0].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[1].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[2].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[0].set_ylabel(\"Depth [m]\")\n",
    "axs[0].set_xlabel(\"Potential density [kg m$^{-3}$]\")\n",
    "axs[1].set_xlabel(r\"$\\epsilon$ [W kg$^{-1}$]\")\n",
    "axs[2].set_xlabel(r\"$N^2$ [rad$^2$ s$^{-2}$]\")\n",
    "\n",
    "fig.align_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The different methods can alter the dissipation rate estimate by almost an order of magnitude. The `endpt` and `teos` methods produce such similar results that they are indistinguishable. The `teosp1` method tends to produce higher $N^2$ estimates, because there is often a jump in stratification adjacent to an overturning patch. The `bulk` method tends to give the smallest $N^2$ estimates. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thorpe scales from temperature alone\n",
    "\n",
    "In some cases it may be desirable to use temperature, rather than density, to detect unstable patches. There are two different ways of doing this:\n",
    "\n",
    "* Provide both temperature and salinity profiles but set `overturns_from_t = True`. In this case, overturns will be diagnosed from the conservative temperature profile or temperature profile, depending on the equation of state chosen. Buoyancy frequency will be estimated from the equation of state or potential density as specified by `N2_method`.\n",
    "\n",
    "* Provide as arguments to `eps_overturn` a profile of temperature and a constant salinity. Pseudo potential density will be calculated assuming a constant salinity profile. Overturns will be diagnosed from the pseudo potential density profile. Depending on the `N2_method` chosen, buoyancy frequency may be estimated from pseudo potential density or the equation of state at constant salinity. As such, buoyancy frequency values will be incorrect in the ocean but perhaps useful for estimating $\\epsilon$. It is probably a good idea to choose a plausible salinity value, such as 0 psu for a fresh water environment or 35 psu for an ocean environment.\n",
    "\n",
    "Note that `dnoise` should describe the temperature noise level if `overturns_from_t = True` and the density noise level at constant salinity if salinity is constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_ip = False\n",
    "N2_method = \"teos\"\n",
    "\n",
    "eps_t, N2_t = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    35,\n",
    "    lon,\n",
    "    lat,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    ")\n",
    "\n",
    "eps_ts, N2_ts = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    dnoise=2e-4,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    "    overturns_from_t=True,\n",
    ")\n",
    "\n",
    "eps, N2 = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    dnoise,\n",
    "    use_ip=use_ip,\n",
    "    N2_method=N2_method,\n",
    ")\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, sharey=True, figsize=(6, 6))\n",
    "axs[0].plot(eps_ts[cut], zc, label=\"CT\")\n",
    "axs[0].plot(eps_t[cut], zc, label=\"pseudo potential density\")\n",
    "axs[0].plot(eps[cut], zc, label=\"potential density (normal case)\")\n",
    "axs[1].plot(N2_ts[cut], zc)\n",
    "axs[1].plot(N2_t[cut], zc)\n",
    "axs[1].plot(N2[cut], zc)\n",
    "axs[0].invert_yaxis()\n",
    "axs[0].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[0].set_ylabel(\"Depth [m]\")\n",
    "axs[0].set_xlabel(r\"$\\epsilon$ [W kg$^{-1}$]\")\n",
    "axs[1].set_xlabel(r\"$N^2$ [rad$^2$ s$^{-2}$]\")\n",
    "\n",
    "fig.align_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The smaller noise parameter in the `overturns_from_t` case means that more overturns are found, however, $\\epsilon$ tends to be a bit lower than the normal case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear equation of state\n",
    "\n",
    "By default, Thorpe scales are estimated from density using the nonlinear [TEOS-10](https://www.teos-10.org/) equation of state. It is also possible to use a linear equation of state, e.g.\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho = \\rho_0 ( 1 - \\alpha (T - T_0) + \\beta (S - S_0)),\n",
    "\\end{equation}\n",
    "\n",
    "which may be useful for numerical modelling or theoretical applications. \n",
    "\n",
    "You may tune the linear parameters as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_EOS_params = dict(rho0=1045.4, a=2e-4, b=7e-4, SP0=35, t0=15)\n",
    "eps_linear, N2_linear, diag_linear = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    use_ip=False,\n",
    "    EOS=\"linear\",\n",
    "    linear_EOS_params=linear_EOS_params,\n",
    "    N2_method=\"bulk\",\n",
    "    return_diagnostics=True,\n",
    ")\n",
    "\n",
    "eps_gsw, N2_gsw, diag_gsw = mx.overturn.eps_overturn(\n",
    "    depth,\n",
    "    t,\n",
    "    SP,\n",
    "    lon,\n",
    "    lat,\n",
    "    use_ip=False,\n",
    "    EOS=\"gsw\",\n",
    "    N2_method=\"bulk\",\n",
    "    return_diagnostics=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot only in the depth range:\n",
    "cut = (depth > 4225) & (depth < 4800)\n",
    "zc = depth[cut]\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, sharey=True, figsize=(9, 6))\n",
    "axs[0].plot(diag_linear[\"dens\"][cut] - 1000, zc, label=\"linear\")\n",
    "# axs[0].plot(diag_linear[\"dens_sorted\"][cut]-1000, zc, label=\"sorted density\")\n",
    "axs[0].plot(diag_gsw[\"dens\"][cut] - 1000, zc, label=\"gsw\")\n",
    "# axs[0].plot(diag_gsw[\"dens_sorted\"][cut]-1000, zc, label=\"sorted density\")\n",
    "axs[1].plot(eps_linear[cut], zc, label=\"linear\")\n",
    "axs[1].plot(eps_gsw[cut], zc, label=\"gsw\")\n",
    "axs[2].plot(N2_linear[cut], zc, label=\"linear\")\n",
    "axs[2].plot(N2_gsw[cut], zc, label=\"gsw\")\n",
    "axs[0].invert_yaxis()\n",
    "axs[0].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[1].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[2].legend(loc=\"lower left\", bbox_to_anchor=(0, 1))\n",
    "axs[0].set_ylabel(\"Depth [m]\")\n",
    "axs[0].set_xlabel(\"Potential density [kg m$^{-3}$]\")\n",
    "axs[1].set_xlabel(r\"$\\epsilon$ [W kg$^{-1}$]\")\n",
    "axs[2].set_xlabel(r\"$N^2$ [rad$^2$ s$^{-2}$]\")\n",
    "\n",
    "fig.align_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how density appears to decrease with depth and the deepest overturn is larger in the linear case. Consider using potential temperature rather than _in situ_ temperature for deep observations when using the linear equation of state, otherwise there is a risk of finding spurious overturns. "
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "notebook_metadata_filter": "-jupytext.text_representation.jupytext_version"
  },
  "kernelspec": {
   "display_name": "Python [conda env:mixsea]",
   "language": "python",
   "name": "conda-env-mixsea-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}