Shear/Strain Parameterization

Note

This section is available as jupyter notebook in the docs/ directory or at https://github.com/modscripps/mixsea/tree/master/docs.

Import mixsea and a few other modules:

import mixsea as mx
import numpy as np
import matplotlib.pyplot as plt

Read example data

ctd = mx.helpers.read_ctd_testfile()
ladcp = mx.helpers.read_ladcp_testfile()

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 for using nan_shearstrain with input data that contains NaNs.

def nonan(data):
    notnan = [np.isfinite(v) for k, v in data.items()]
    notnan = np.vstack(notnan)
    notnan = np.all(notnan, axis=0)
    return notnan


notnan = nonan(ctd)
depth = ctd["depth"][notnan]
t = ctd["t"][notnan]
SP = ctd["SP"][notnan]
lon = ctd["lon"][0]
lat = ctd["lat"][0]

notnan = nonan(ladcp)
u = ladcp["u"][notnan]
v = ladcp["v"][notnan]
uz = ladcp["uz"][notnan]
vz = ladcp["vz"][notnan]
depth_sh = ladcp["depth"][notnan]

A quick overview plot of the data:

fig, ax = plt.subplots(
    nrows=1, ncols=3, figsize=(9, 4), constrained_layout=True, sharey=True
)
ax[0].plot(t, depth)
ax[0].set(ylabel="depth [m]", xlabel="temperature [°C]")
ax[1].plot(u, depth_sh, label="u")
ax[1].plot(v, depth_sh, label="v")
ax[1].set(xlabel="velocity [m/s]")
ax[1].legend()
ax[2].plot(uz, depth_sh, label=r"u$_{z}$")
ax[2].plot(vz, depth_sh, label=r"v$_{z}$")
ax[2].set(xlabel="shear [1/s]")
ax[2].legend()
ax[0].invert_yaxis()

Stratification Smoothing Techniques

shearstrain allows to select between two different techniques for obtaining a smooth profile of \(N^2\).

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).

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).

Set up shear/strain input parameters:

shst_params = dict()
# Center points of depth windows. Windows are half overlapping, i.e.
# their size (300m) is double the spacing here (150m).
window_size = 300
dz = window_size / 2
print("window size {} m, window spacing {} m".format(window_size, dz))
shst_params["depth_bin"] = np.linspace(dz, dz * 80, num=80)
shst_params["window_size"] = window_size
# Wavenumber vector.
shst_params["m"] = np.arange(
    2 * np.pi / window_size, 2 * np.pi / 10, 2 * np.pi / window_size
)
# Wavenumber indices for integration. Shear is integrated from 300m to
# 100m scales. Strain is integrated from 150m to 30m.
shst_params["m_include_sh"] = list(range(3))
shst_params["m_include_st"] = list(range(1, 12))

shst_params["ladcp_is_shear"] = True
shst_params["return_diagnostics"] = True

window size 300 m, window spacing 150.0 m

Run the shear/strain parameterization for both \(N^2\) techniques:

eps_shst, krho_shst, diag = mx.shearstrain.shearstrain(
    depth, t, SP, lon, lat, uz, vz, depth_sh, smooth="PF", **shst_params
)

eps_shst2, krho_shst2, diag2 = mx.shearstrain.shearstrain(
    depth, t, SP, lon, lat, uz, vz, depth_sh, smooth="AL", **shst_params
)

depth_bin = diag["depth_bin"]
fig, ax = plt.subplots(
    nrows=1, ncols=2, figsize=(6, 4), constrained_layout=True, sharey=True
)
opts = dict(drawstyle="steps-mid")
# shear/strain
ax[0].plot(eps_shst, depth_bin, label="shear/strain PF", **opts)
ax[0].plot(eps_shst2, depth_bin, label="shear/strain AL", **opts)
ax[0].legend()
ax[0].set(
    xscale="log", xlabel=r"$\epsilon$ [W/kg]", ylabel="depth [m]",
)
# strain only
eps_st = diag["eps_st"]
eps_st2 = diag2["eps_st"]
ax[1].plot(eps_st, depth_bin, label="strain only PF", **opts)
ax[1].plot(eps_st2, depth_bin, label="strain only AL", **opts)
ax[1].legend()
ax[1].set(
    xscale="log", xlabel=r"$\epsilon$ [W/kg]",
)
ax[0].invert_yaxis()

Missing Data

Here we pass the data as is, containing NaNs, to nan_shearstrain.

eps_shst3, krho_shst3, diag3 = mx.shearstrain.nan_shearstrain(
    ctd["depth"],
    ctd["t"],
    ctd["SP"],
    ctd["lon"],
    ctd["lat"],
    ladcp["uz"],
    ladcp["vz"],
    ladcp["depth"],
    smooth="PF",
    **shst_params
)

Compare with result of shearstrain from above - the results should be the same:

depth_bin = diag["depth_bin"]
depth_bin3 = diag3["depth_bin"]
fig, ax = plt.subplots(
    nrows=1, ncols=2, figsize=(6, 4), constrained_layout=True, sharey=True
)
opts = dict(drawstyle="steps-mid")
# shear/strain
ax[0].plot(eps_shst, depth_bin, label="shear/strain PF", **opts)
ax[0].plot(eps_shst3, depth_bin3, label="shear/strain PF", linestyle="--", **opts)
ax[0].legend()
ax[0].set(
    xscale="log", xlabel=r"$\epsilon$ [W/kg]", ylabel="depth [m]",
)
# strain only
eps_st = diag["eps_st"]
eps_st3 = diag3["eps_st"]
ax[1].plot(eps_st, depth_bin, label="strain only PF", **opts)
ax[1].plot(eps_st3, depth_bin3, label="strain only PF", linestyle="--", **opts)
ax[1].legend()
ax[1].set(
    xscale="log", xlabel=r"$\epsilon$ [W/kg]",
)
ax[0].invert_yaxis()