docs/shearstrain.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "52ba291eebef",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "source": [
    "# Shear/Strain Parameterization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "f2fd47f80438b",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "This section is available as jupyter notebook in the `docs/` directory or at https://github.com/modscripps/mixsea/tree/main/docs.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "503fd9cd7e3d18",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "source": [
    "Import `mixsea` and a few other modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "5095fbea004c38",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "outputs": [],
   "source": [
    "import mixsea as mx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import gsw\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Basic idea**\n",
    "\n",
    "(Unfortunately shear-strain parameterizations for internal wave-driven mixing are not basic at all; this is as clear and concise a summary as possible.)\n",
    "\n",
    "At the root of the shear-strain parameterization is the assumption of a scale separation between the production of turbulent kinetic energy $\\mathcal{P}$ (extracted from a mean shear) at low wave numbers and its destruction by both viscous dissipation $\\epsilon$ and vertical buoyancy flux $\\overline{w'b'}$ at high wave numbers. In stationary turbulence, non-linear wave-wave interactions at some intermediate (\"fine\"-)scale $m_{c}$ must be responsible for this downscale energy flux $\\mathcal{T}_{\\text{IW}}$, i.e.\n",
    "\\begin{equation}\n",
    "-\\mathcal{T}_{\\text{IW}}(m_{c}) =  \\mathcal{P} = - \\epsilon + \\overline{w'b'}\n",
    "\\end{equation}\n",
    "\n",
    "Shear-strain parameterizations get their popular name based on the functional form $\\mathcal{G}$ of the downscale transfer function that arises from either heuristic (<cite data-cite=\"Polzin2004\">Polzin 2004</cite>) or ray-tracing (<cite data-cite=\"Henyey1986\">Henyey et al. 1986</cite>) methods:\n",
    "\\begin{equation}\n",
    "\\mathcal{T}_{\\text{IW}} = \\mathcal{G}(E_{\\text{KE}}, E_{\\text{PE}}; \\overline{N}, f),\n",
    "\\end{equation}\n",
    "where:\n",
    "\n",
    "- $E_{\\text{KE}}$ is the *kinetic energy spectral density*, often expressed as the shear variance $\\langle S^{2} \\rangle = \\int_{m_{\\text{low}}}^{\\min\\{m_{\\text{high}}, m_{c}\\}} m^{2}E_{\\text{KE}} \\; \\text{d} m$ integrated over a low wavenumber band [$m_{\\text{low}}$, $m_{\\text{high}}$], which is possibly cutoff early at $m_{c}$ as soon the cumulative shear variance increases above $1/6$ the amount required for wave breaking by shear instability $\\big(\\text{i.e. }R_{f} = \\langle N^{2} \\rangle / \\langle S^{2} \\rangle = 1/4 \\big)$;\n",
    "    \n",
    "    \n",
    "- $E_{\\text{PE}}$ is the *potential energy spectral density*, often expressed in terms of a strain variance $\\langle \\xi_{z}^{2} \\rangle$, where $\\xi$ is the displacement of a density surface, or in terms of the ratio  of shear to strain variance, $R_{\\omega} = E_{\\text{KE}}/E_{\\text{PE}}$, which also doubles as a characteristic internal wave aspect ratio (strain variance also has a high cutoff wavenumber);\n",
    "    \n",
    "    \n",
    "- $\\overline{N}$ is the mean stratification; and \n",
    "    \n",
    "    \n",
    "- $f$ is the Coriolis frequency.\n",
    "    \n",
    "    \n",
    "Strain-only parameterizations assume a characteristic constant shear to strain ratio and thus replace the shear variance with more readily available strain variance measurements, $\\langle S^{2} \\rangle \\simeq \\langle \\xi_{z}^{2} \\rangle R_{\\omega}$.\n",
    "\n",
    "In practice, such parameterizations are conventionally formulated as scalings of the <cite data-cite=\"Garrett1975\">Garrett-Munk (1975)</cite> empirical internal wave energy spectrum. The exact forms of the parameterization used are described in the underlying functions' docstrings– see the Methods section of <cite data-cite=\"Gregg2003\">Gregg (2003)</cite> for a concise example or <cite data-cite=\"Kunze2006\">Kunze (2006)</cite> for a detailed exposition of the default method implemented in `mixsea`.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read example data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "06bb7b7703f98",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "outputs": [],
   "source": [
    "ctd = mx.helpers.read_ctd_testfile()\n",
    "ladcp = mx.helpers.read_ladcp_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. See the section [Missing Data](#Missing-Data) for using `nan_shearstrain` with input data that contains NaNs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonan(data):\n",
    "    notnan = [np.isfinite(v) for k, v in data.items()]\n",
    "    notnan = np.vstack(notnan)\n",
    "    notnan = np.all(notnan, axis=0)\n",
    "    return notnan\n",
    "\n",
    "\n",
    "notnan = nonan(ctd)\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]\n",
    "\n",
    "notnan = nonan(ladcp)\n",
    "u = ladcp[\"u\"][notnan]\n",
    "v = ladcp[\"v\"][notnan]\n",
    "uz = ladcp[\"uz\"][notnan]\n",
    "vz = ladcp[\"vz\"][notnan]\n",
    "depth_sh = ladcp[\"depth\"][notnan]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "7647faceff63",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "source": [
    "A quick overview plot of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "f5fe8d01d3c",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(\n",
    "    nrows=1, ncols=3, figsize=(9, 4), constrained_layout=True, sharey=True\n",
    ")\n",
    "ax[0].plot(t, depth)\n",
    "ax[0].set(ylabel=\"depth [m]\", xlabel=\"temperature [°C]\")\n",
    "ax[1].plot(u, depth_sh, label=\"u\")\n",
    "ax[1].plot(v, depth_sh, label=\"v\")\n",
    "ax[1].set(xlabel=\"velocity [m/s]\")\n",
    "ax[1].legend()\n",
    "ax[2].plot(uz, depth_sh, label=r\"u$_{z}$\")\n",
    "ax[2].plot(vz, depth_sh, label=r\"v$_{z}$\")\n",
    "ax[2].set(xlabel=\"shear [1/s]\")\n",
    "ax[2].legend()\n",
    "ax[0].invert_yaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example of the shear/strain internal wave parameterization\n",
    "\n",
    "We begin by choosing our parameters. The key parameter choices are:\n",
    "\n",
    "- `window_size`: the vertical extent of the spectral windows used to estimate shear/strain spectra;\n",
    "- `min_size`: the shortest vertical wavelength resolved by the data (should be a factor of `window_size`); doesn't matter much in practice because we separately constrain the spectral windows with the two following parameters.\n",
    "- `mi_sh`: index limits for integrating shear variance (excluding upper limit, following Python convention). Lower bound is generally just set to 0 so that the gravest mode included is the `window_size`; upper bound is generally set to 100-150 m.\n",
    "- `mi_st`: index limits for integrating strain variance (excluding upper limit, following Python convention). Lower bound is generally set to 100-150 m and upper bound set to 10-20 m. In practice, the integration is often cut-off well before 10 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shst_params = dict()\n",
    "\n",
    "# Center points of depth windows. Windows are half overlapping, i.e.\n",
    "# their size (300m) is double the spacing here (150m).\n",
    "window_size = 300.0\n",
    "min_size = 10.0\n",
    "dz = window_size / 2\n",
    "shst_params[\"depth_bin\"] = np.arange(dz, 10000.0, dz)\n",
    "shst_params[\"window_size\"] = window_size\n",
    "\n",
    "# Set up wavenumber vector.\n",
    "shst_params[\"m\"] = np.arange(\n",
    "    2 * np.pi / window_size, 2 * np.pi / min_size, 2 * np.pi / window_size\n",
    ")\n",
    "\n",
    "# Set up limits for shear and strain variance integrations\n",
    "mi_sh = np.array([0, 3])\n",
    "mii_sh = np.array(range(*mi_sh))\n",
    "mi_st = np.array([2, 20])\n",
    "mii_st = np.array(range(*mi_st))\n",
    "\n",
    "shst_params[\"m_include_sh\"] = mii_sh\n",
    "shst_params[\"m_include_st\"] = mii_st\n",
    "\n",
    "# Convert indices to more intuitive length scales\n",
    "m_sh = 2 * np.pi / shst_params[\"m\"][[mi_sh[0], mi_sh[1] - 1]]\n",
    "m_st = 2 * np.pi / shst_params[\"m\"][[mi_st[0], mi_st[1] - 1]]\n",
    "print(\n",
    "    f\"Wavenumber indices for integration:\\n\"\n",
    "    f\"- Shear is integrated from {round(m_sh[0])}m to {round(m_sh[1])}m scales.\\n\"\n",
    "    f\"- Strain is integrated from {round(m_st[0])}m to {round(m_st[1])}m.\"\n",
    ")\n",
    "\n",
    "shst_params[\"ladcp_is_shear\"] = True\n",
    "shst_params[\"return_diagnostics\"] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now pass our input data and parameter dictionary to the `shearstrain` function and retrieve the dissipation rate of turbulent kinetic energy $\\epsilon$, the diapycnal diffusivity $\\kappa_{\\rho}$, and an optional dictionary of detailed diagnostics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps, krho, diag = mx.shearstrain.shearstrain(\n",
    "    depth, t, SP, lon, lat, uz, vz, depth_sh, **shst_params\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the vertical profiles of the estimated mixing rates due to internal waves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_bin = diag[\"depth_bin\"]\n",
    "depth_bin_edges = np.concatenate(\n",
    "    (\n",
    "        [np.min(depth_bin) - dz / 2],\n",
    "        0.5 * (depth_bin[1:] + depth_bin[:-1]),\n",
    "        [np.max(depth_bin) + dz / 2],\n",
    "    )\n",
    ")\n",
    "fig, ax = plt.subplots(\n",
    "    nrows=1, ncols=2, figsize=(6, 4), sharey=True, constrained_layout=True\n",
    ")\n",
    "ax[0].stairs(eps, depth_bin_edges, orientation=\"horizontal\", baseline=None, lw=1.5)\n",
    "ax[0].set(xscale=\"log\", xlabel=r\"$\\epsilon$ [W/kg]\", ylabel=\"depth [m]\")\n",
    "ax[0].invert_yaxis()\n",
    "\n",
    "ax[1].stairs(krho, depth_bin_edges, orientation=\"horizontal\", baseline=None, lw=1.5)\n",
    "ax[1].set(xscale=\"log\", xlabel=r\"$\\kappa$ [m$^{2}$/s]\")\n",
    "\n",
    "k = -2\n",
    "ax[0].fill_betweenx(\n",
    "    np.array([-dz, dz]) + depth_bin[k],\n",
    "    [0, 0],\n",
    "    [eps[k], eps[k]],\n",
    "    facecolor=\"C3\",\n",
    "    alpha=0.2,\n",
    ")\n",
    "ax[1].fill_betweenx(\n",
    "    np.array([-dz, dz]) + depth_bin[k],\n",
    "    [0, 0],\n",
    "    [krho[k], krho[k]],\n",
    "    facecolor=\"C3\",\n",
    "    alpha=0.2,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Deep dive for a single high-dissipation window near the bottom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's zoom in on the second to last window (shaded in red above), which exhibits abnormally strong strain variance and thus corresponds to a large dissipation rate, to get a better understanding of what the shear/strain variance and strain variance algorithms are doing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot()\n",
    "plt.loglog(diag[\"m\"][mii_sh], diag[\"P_shear\"][k][mii_sh])\n",
    "plt.loglog(diag[\"m\"], diag[\"P_shear\"][k], \"C0--\", alpha=0.25)\n",
    "S_shear_gm, P_shear_gm = mx.shearstrain.gm_shear_variance(\n",
    "    diag[\"m\"], mii_sh, diag[\"Nmseg\"][k]\n",
    ")\n",
    "plt.loglog(diag[\"m\"], P_shear_gm, label=\"GM76 spectrum (minus rolloff)\")\n",
    "ylims = plt.gca().get_ylim()\n",
    "plt.fill_between(\n",
    "    diag[\"m\"][[mi_sh[0], mi_sh[1] - 1]],\n",
    "    [ylims[0], ylims[0]],\n",
    "    [ylims[1], ylims[1]],\n",
    "    alpha=0.1,\n",
    ")\n",
    "plt.ylim(ylims)\n",
    "plt.xlim(diag[\"m\"][0], diag[\"m\"][-1])\n",
    "plt.ylabel(r\"$\\langle S^{2} \\rangle$\", fontsize=13)\n",
    "plt.xlabel(r\"wavenumber $m$ [rad/m]\", fontsize=13)\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "def same(x):\n",
    "    return x\n",
    "\n",
    "\n",
    "secax = ax.secondary_xaxis(\"top\", functions=(same, same))\n",
    "secax.set_xticks([diag[\"m\"][mii_sh[0]], diag[\"m\"][mii_sh[-1]]], minor=False)\n",
    "secax.set_xticklabels([r\"$m_{low}$\", r\"$m_{high}$\"])\n",
    "\n",
    "var_ratio = diag[\"Int_sh\"][k] / diag[\"Int_sh_gm\"][k]\n",
    "print(\n",
    "    f\"Shear variance levels are {round(var_ratio,1)} times higher than\\n\"\n",
    "    f\"the GM spectrum over the low wavenumber band.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot()\n",
    "plt.loglog(diag[\"m\"][mii_st], diag[\"P_strain\"][k][mii_st])\n",
    "plt.loglog(diag[\"m\"], diag[\"P_strain\"][k], \"C0--\", alpha=0.25)\n",
    "S_strain_gm, P_strain_gm = mx.shearstrain.gm_strain_variance(\n",
    "    diag[\"m\"], mii_st, diag[\"Nmseg\"][k]\n",
    ")\n",
    "plt.loglog(diag[\"m\"], P_strain_gm, label=\"GM76 spectrum\")\n",
    "mii_st_cutoff = mii_st[diag[\"m\"][mii_st] <= diag[\"Mmax_st\"][k]]\n",
    "S_strain_gm, P_strain_gm = mx.shearstrain.gm_strain_variance(\n",
    "    diag[\"m\"], mii_st_cutoff, diag[\"Nmseg\"][k]\n",
    ")\n",
    "ylims = plt.gca().get_ylim()\n",
    "plt.fill_between(\n",
    "    [diag[\"m\"][mii_st_cutoff[0]], diag[\"Mmax_st\"][k]],\n",
    "    [ylims[0], ylims[0]],\n",
    "    [ylims[1], ylims[1]],\n",
    "    alpha=0.1,\n",
    ")\n",
    "plt.ylim(ylims)\n",
    "plt.xlim(diag[\"m\"][0], diag[\"m\"][-1])\n",
    "plt.ylabel(r\"$\\langle \\xi_{z}^{2} \\rangle$\", fontsize=13)\n",
    "plt.xlabel(r\"wavenumber $m$ [rad/m]\", fontsize=13)\n",
    "plt.legend()\n",
    "\n",
    "secax = ax.secondary_xaxis(\"top\", functions=(same, same))\n",
    "secax.set_xticks(\n",
    "    [diag[\"m\"][mii_st[0]], diag[\"m\"][mii_st[-1]], diag[\"Mmax_st\"][k]], minor=False\n",
    ")\n",
    "secax.set_xticklabels([r\"$m_{low}$\", r\"$m_{high}$\", r\"$m_{c}$\"])\n",
    "\n",
    "print(\n",
    "    rf\"Strain variance levels are {round(diag['Int_st'][k]/S_strain_gm,1)} times \"\n",
    "    f\"higher than\\nthe GM spectrum over the low wavenumber band.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, since the strain variance is much larger than GM levels while the shear variance is only slightly large than GM levels, the shear/strain ratio drops well below the GM level of $R_{\\omega}$=3. Since shear/strain ratios of less than $1$ are unphysical for internal waves (for example the aspect ratio correction terms contain a factor $\\sqrt{R_{\\omega} - 1}$), we set any $R_{\\omega} \\leq 1$ to $R_{\\omega}=1.01$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"Raw shear/strain ratio: {round(diag['Rwtot'][k], 2)}, \"\n",
    "    f\"corrected ratio: {diag['Rwcor'][k]}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Decomposing the terms of the parameterization\n",
    "\n",
    "The dissipation rate for the shear/strain parameterization is given as:\n",
    "\\begin{equation}\n",
    "\\epsilon_{\\text{fine}} = \\epsilon_{0} \\; \\frac{\\overline{N}^{2}}{N_{0}^{2}} \\; \\frac{\\langle S^{2} \\rangle}{\\langle S^{2}_{\\text{GM}} \\rangle} \\; h_{1}(R_{\\omega}) \\; L(f, \\overline{N})\n",
    "\\end{equation}\n",
    "\n",
    "It is useful to plot their vertical profiles to get a sense of how their various dependencies combine to produce a typical dissipation rate profile, where the grey dashed lines denote the reference values for the GM spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N0 = 5.24e-3\n",
    "f = np.absolute(gsw.f(ctd[\"lat\"][0]))\n",
    "\n",
    "fig, ax = plt.subplots(\n",
    "    nrows=1, ncols=5, figsize=(9, 4), sharey=True, constrained_layout=True\n",
    ")\n",
    "ax[0].stairs(\n",
    "    eps,\n",
    "    depth_bin_edges,\n",
    "    label=r\"$\\epsilon$\",\n",
    "    orientation=\"horizontal\",\n",
    "    baseline=None,\n",
    "    lw=1.5,\n",
    ")\n",
    "ax[0].set(xscale=\"log\", xlabel=r\"$\\epsilon$ [W/kg]\", ylabel=\"depth [m]\")\n",
    "ax[0].invert_yaxis()\n",
    "\n",
    "ax[1].axvline(1.0, color=\"k\", alpha=0.3, ls=\"--\")\n",
    "ax[1].stairs(\n",
    "    diag[\"Nmseg\"] ** 2 / N0**2,\n",
    "    depth_bin_edges,\n",
    "    orientation=\"horizontal\",\n",
    "    baseline=None,\n",
    "    lw=1.5,\n",
    ")\n",
    "ax[1].set(xscale=\"log\", xlabel=r\"$\\overline{N}^{2}/N_{0}^{2}$\")\n",
    "\n",
    "ax[2].axvline(1.0, color=\"k\", alpha=0.3, ls=\"--\")\n",
    "ax[2].stairs(\n",
    "    diag[\"Int_sh\"] ** 2 / diag[\"Int_sh_gm\"] ** 2,\n",
    "    depth_bin_edges,\n",
    "    orientation=\"horizontal\",\n",
    "    baseline=None,\n",
    "    lw=1.5,\n",
    ")\n",
    "ax[2].set(\n",
    "    xscale=\"log\",\n",
    "    xlabel=r\"$\\langle S^{2} \\rangle^{2} / \\langle S^{2}_{GM} \\rangle^{2}$ [W/kg]\",\n",
    ")\n",
    "\n",
    "ax[3].axvline(\n",
    "    mx.shearstrain.aspect_ratio_correction_shst(3.0), color=\"k\", alpha=0.3, ls=\"--\"\n",
    ")\n",
    "ax[3].stairs(\n",
    "    mx.shearstrain.aspect_ratio_correction_shst(diag[\"Rwcor\"]),\n",
    "    depth_bin_edges,\n",
    "    orientation=\"horizontal\",\n",
    "    baseline=None,\n",
    "    lw=1.5,\n",
    ")\n",
    "ax[3].set(xscale=\"log\", xlabel=r\"$h_{1}(R_{\\omega})$\")\n",
    "\n",
    "ax[4].axvline(1.0, color=\"k\", alpha=0.3, ls=\"--\")\n",
    "ax[4].stairs(\n",
    "    mx.shearstrain.latitude_correction(f, diag[\"Nmseg\"]),\n",
    "    depth_bin_edges,\n",
    "    orientation=\"horizontal\",\n",
    "    baseline=None,\n",
    "    lw=1.5,\n",
    ")\n",
    "ax[4].set(xlabel=r\"$L(f, N)$\")\n",
    "\n",
    "eps0 = 7.8e-10\n",
    "eps_verify = (\n",
    "    eps0\n",
    "    * diag[\"Nmseg\"] ** 2\n",
    "    / N0**2\n",
    "    * diag[\"Int_sh\"] ** 2\n",
    "    / diag[\"Int_sh_gm\"] ** 2\n",
    "    * mx.shearstrain.aspect_ratio_correction_shst(diag[\"Rwcor\"])\n",
    "    * mx.shearstrain.latitude_correction(f, diag[\"Nmseg\"])\n",
    ")\n",
    "ax[0].stairs(\n",
    "    eps_verify,\n",
    "    depth_bin_edges,\n",
    "    label=r\"$\\epsilon_{verify}$\",\n",
    "    orientation=\"horizontal\",\n",
    "    baseline=None,\n",
    "    lw=1.2,\n",
    "    ls=\"--\",\n",
    "    alpha=1.0,\n",
    ")\n",
    "ax[0].legend()\n",
    "\n",
    "print(f\"Verification successful: {np.isclose(eps_verify, eps).all()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "640def2242cea8",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "source": [
    "## Stratification Smoothing Techniques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`shearstrain` allows to select between two different techniques for obtaining a smooth profile of $N^2$. Adiabatic leveling is the default (used in the example above).\n",
    "\n",
    "**Adiabatic Leveling** finds locally stable density gradients to calculate a smooth $N^2$ profile following Bray & Fofonoff (1981), see e.g. Polzin et al. (2014), Waterman et al. (2013, 2014).\n",
    "\n",
    "**Polynomial Fitting** applies a second order polynomial fit to $N^2$ for each segment to obtain a smooth profile, compare e.g. Kunze et al. (2006), Whalen et al. (2012, 2015)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "c41f8f269f8438",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "source": [
    "Run the shear/strain parameterization for both the techniques:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "83b2ab34a7137",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "outputs": [],
   "source": [
    "eps_PF, krho_PF, diag_PF = mx.shearstrain.shearstrain(\n",
    "    depth, t, SP, lon, lat, uz, vz, depth_sh, smooth=\"PF\", **shst_params\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare results from the two techniques:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "ce67995fee7888",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(\n",
    "    nrows=1, ncols=2, figsize=(6, 4), constrained_layout=True, sharey=True\n",
    ")\n",
    "opts = dict(orientation=\"horizontal\", baseline=None)\n",
    "\n",
    "# shear/strain\n",
    "ax[0].stairs(eps, depth_bin_edges, label=\"AL\", **opts)\n",
    "ax[0].stairs(eps_PF, depth_bin_edges, label=\"PF\", **opts)\n",
    "ax[0].legend(loc=\"upper left\")\n",
    "ax[0].set(\n",
    "    xscale=\"log\", xlabel=r\"$\\epsilon$ [W/kg]\", ylabel=\"depth [m]\", title=\"shear/strain\"\n",
    ")\n",
    "ax[0].invert_yaxis()\n",
    "\n",
    "# strain only\n",
    "ax[1].stairs(diag[\"eps_st\"], depth_bin_edges, label=\"AL\", **opts)\n",
    "ax[1].stairs(diag_PF[\"eps_st\"], depth_bin_edges, label=\"PF\", **opts)\n",
    "ax[1].set(xscale=\"log\", xlabel=r\"$\\epsilon$ [W/kg]\", title=\"strain only\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results differ mostly for the upper ocean and the bottom layer. Note that we do not filter out parts of the profile where the shear/strain parameterization may not apply as underlying assumptions are violated, i.e. near the surface mixed layer and in the lowly stratified layer of the abyssal flow through the Samoan Passage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Missing-Data'></a>\n",
    "\n",
    "## Missing Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we pass the data as is, containing NaNs, to `nan_shearstrain`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps_nan, krho_nan, diag_nan = mx.shearstrain.nan_shearstrain(\n",
    "    ctd[\"depth\"],\n",
    "    ctd[\"t\"],\n",
    "    ctd[\"SP\"],\n",
    "    ctd[\"lon\"],\n",
    "    ctd[\"lat\"],\n",
    "    ladcp[\"uz\"],\n",
    "    ladcp[\"vz\"],\n",
    "    ladcp[\"depth\"],\n",
    "    **shst_params\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with result of `shearstrain` from above - the results should be the same. Here we apply only the polynomial fit method for obtaining strain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "janus": {
     "all_versions_showing": false,
     "cell_hidden": false,
     "current_version": 0,
     "id": "ce67995fee7888",
     "named_versions": [],
     "output_hidden": false,
     "show_versions": false,
     "source_hidden": false,
     "versions": []
    }
   },
   "outputs": [],
   "source": [
    "depth_bin_nan = diag_nan[\"depth_bin\"]\n",
    "depth_bin_nan_edges = np.concatenate(\n",
    "    (\n",
    "        [np.min(depth_bin_nan) - dz / 2],\n",
    "        0.5 * (depth_bin_nan[1:] + depth_bin_nan[:-1]),\n",
    "        [np.max(depth_bin_nan) - dz / 2],\n",
    "    )\n",
    ")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(\n",
    "    nrows=1, ncols=2, figsize=(6, 4), constrained_layout=True, sharey=True\n",
    ")\n",
    "\n",
    "# shear/strain\n",
    "ax[0].stairs(eps, depth_bin_edges, label=\"without NaNs\", **opts)\n",
    "ax[0].stairs(eps_nan, depth_bin_nan_edges, label=\"with NaNs\", linestyle=\"--\", **opts)\n",
    "ax[0].legend()\n",
    "ax[0].set(\n",
    "    xscale=\"log\", xlabel=r\"$\\epsilon$ [W/kg]\", ylabel=\"depth [m]\", title=\"shear/strain\"\n",
    ")\n",
    "ax[0].invert_yaxis()\n",
    "\n",
    "# strain only\n",
    "ax[1].stairs(diag[\"eps_st\"], depth_bin_edges, label=\"without NaNs\", **opts)\n",
    "ax[1].stairs(\n",
    "    diag_nan[\"eps_st\"], depth_bin_nan_edges, label=\"with NaNs\", linestyle=\"--\", **opts\n",
    ")\n",
    "ax[1].legend()\n",
    "ax[1].set(\n",
    "    xscale=\"log\",\n",
    "    xlabel=r\"$\\epsilon$ [W/kg]\",\n",
    "    title=\"strain only\",\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there may be cases where simply removing all missing data from the input, as is happening in `nan_shearstrain`, will lead to unphysical results! "
   ]
  }
 ],
 "metadata": {
  "janus": {
   "filepaths": [],
   "janus_markers": [],
   "track_history": true
  },
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-"
  },
  "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
}