Finescale Mixing Parametrizations: Development

Finescale Mixing Parametrizations: Development#

From Whalen et al (2012) Supplementary

QC:

All Argo float profiles with an ‘A’ quality rating (all real-time quality control tests passed) between June 2006 and December 2011 were selected.
The mixed layer was removed before implementing the parameterization since the mixed layer appears at a region of the high strain variance, and the parameterization treats all strain variance as internal waves, leading to an inaccurate estimate of the mixing. The same problem exists for other areas of low stratification near the surface (for example, mode water).
The variable temperature criterion (de Boyer Montegut et al., 2004) was therefore applied twice, once to remove the mixed layer, and a second time treating the bottom of the mixed layer as the surface.
Profiles with large spikes from any of the sensors were taken out of the dataset, and small salinity spikes needed to be removed from a subset of the profiles.

Steps:

Each profile containing 15 m resolution or better was divided into 200 m half-overlapping segments.
Segments were discarded if N² < 10−9 s−2 or if N² varied by more than 6 x 10−4 s−2, to remove sharp strain gradients not associated with internal waves.
We consider the vertical resolution of the segments since this could have an effect on the calculated dissipation rates. The majority of the segments had either 2 m or 10 m resolution, with similar mean ! (the mean 10 m resolution estimate above 500 m was only 10% larger), while the standard deviation of dissipation rate for 2 m data was 40% larger.

Whalen et al (2015):

To calculate the strain variance ⟨ξ²⟩, each segment was detrended, windowed using a sin² 10% taper, and spatially Fourier transformed to generate the spectrum Sstr for each segment. The spectrum is corrected for first differencing by dividing by the transfer function sinc² (kzDz/2p), where Dz is the vertical resolution of the segment. Argo data are either point measurements or averages over a depth interval, which is currently not identified in the metadata. Here, we correct for this whenever we know the sampling scheme of the profile by dividing by the same transfer function a second time [a correction not included in Whalen et al. (2012)]. This slightly raises the variances in the bin-averaged cases. The average increase in the dissipation rate for the Atlantic Ocean is a factor of 1.02, with a range of 1 to 2.6.

Strain \begin{equation} ξ_z = \frac{N² - N²_{fit}}{\bar{N²}} \end{equation}

Waterman et al (2014)

Segments were constructed both starting from the surface and the bottom, with near-surface values computed from the vertical segments defined starting from the surface, near-bottom values computed from the vertical segments defined starting from the bottom, and interior values computed from an average of top–down- and bottom–up-defined segments.

%load_ext watermark

import argopy
import cf_xarray as cfxr
import dask
import dcpy
import gsw
import matplotlib as mpl
import matplotlib.pyplot as plt
import mixsea
import numpy as np
import pandas as pd
import scipy
import scipy.integrate
import seawater as sw
from dcpy.finestructure import estimate_turb_segment
from scipy import signal
from scipy.io import loadmat

import xarray as xr

plt.style.use("bmh")
mpl.rcParams["figure.dpi"] = 140
mpl.rcParams["lines.linewidth"] = 1
dirname = "/home/deepak/datasets/finestructure/"

%watermark -iv

dask      : 2021.6.2
cf_xarray : 0.4.1.dev21+gab9dc66
scipy     : 1.5.3
eddydiff  : 0.1
pandas    : 1.2.3
seawater  : 3.3.4
dcpy      : 0.1
mixsea    : 0.1.0
xarray    : 0.17.1.dev3+g48378c4b1
numpy     : 1.20.2
argopy    : 0.1.7
gsw       : 3.4.0
matplotlib: 3.4.1

kunze = dcpy.oceans.read_kunze_2017_finestructure()

Check against Kunze & mixsea#

cruise = "317519930704"
cast = 19
kunze_cruise = kunze.groupby("cruise").get_group(cruise).query(f"drop == {cast}")

ctd = xr.open_mfdataset(f"{dirname}/../cchdo/{cruise}/*ctd.nc").load()
del ctd.attrs["comments"]

profile = ctd.query({"N_PROF": f"cast == {cast}"}).squeeze()
profile.ctd_salinity.attrs["standard_name"] = "sea_water_salinity"
profile["γ"] = dcpy.oceans.neutral_density(profile).reset_coords(drop=True)

# Polzin et al (2015) recommend sorting then calculating N2
initial_pressure = profile.pressure.reset_coords(drop=True)
profile = profile.sortby("γ")
profile["pressure"] = initial_pressure
profile = profile.swap_dims({"N_LEVELS": "pressure"})
profile = profile.isel(pressure=profile.pressure.notnull())

profile["N2"] = (
    "pressure",
    np.append(
        dcpy.eos.bfrq(
            profile.ctd_salinity,
            profile.ctd_temperature,
            profile.pressure,
            dim="pressure",
        )[0]
        # .interp(pressure_mid=profile.pressure.data)
        .data[1:],
        [np.nan, np.nan],
    ),
)

SA = gsw.SA_from_SP(
    profile.ctd_salinity, profile.pressure, profile.longitude, profile.latitude
)
CT = gsw.CT_from_t(SA, profile.ctd_temperature, profile.pressure)
N2gsw, p = gsw.Nsquared(SA, CT, profile.pressure, profile.latitude)
N2gsw = xr.DataArray(N2gsw, dims="pressure", coords={"pressure": p})

# profile["N2gsw"] = ("pressure", N2gsw.interp(pressure=profile.pressure.data).data)
profile["N2gsw"] = ("pressure", np.append(N2gsw, np.nan))

# Density variables σp, σθ, σ3, and γ_n and stratification N2 were then computed
profile["σθ"] = dcpy.eos.pden(
    profile.ctd_salinity, profile.ctd_temperature, profile.pressure, 0
)
profile["σ3"] = dcpy.eos.pden(
    profile.ctd_salinity, profile.ctd_temperature, profile.pressure, 3000
)

profile["N2θ"] = (
    -9.81
    / 1030
    * profile.σθ.cf.differentiate("sea_water_pressure", positive_upward=True)
)
profile["N23"] = (
    -9.81
    / 1030
    * profile.σ3.cf.differentiate("sea_water_pressure", positive_upward=True)
)
profile["N2γ"] = (
    -9.81
    / 1030
    * profile.γ.cf.differentiate("sea_water_pressure", positive_upward=True)
)

profile["pres"] = profile["pressure"]

Break up in to half overlapping segments

Compare to mixsea & Sensitivity studies#

The real difference is because they do some funny interpolating of the spectrum that I do not understand: https://github.com/modscripps/mixsea/blob/2c1a27b5b18d8b28f1af17e6d7dbe0f505ef78d8/mixsea/helpers.py#L169-L178

If I comment out that interp step, spectra agree exactly.

Sensitivity

The problem is that the cutoff wavenumber chosen depends a lot on spectrum; and then that determines how much of the GM spectrum is integrated. Taken together you get a lot of sensitivity just from that :/
Method of choosing the cutoff wavenumber matters too
So does Rω or shear-strain ratio.

I think part of the disagreement is also because the mixsea thing isn’t using the same bins

seg = segmented.isel(pressure=9)

dp = 2
N2 = seg.N2gsw
seg["N2fit"] = (
    seg.N2.dims,
    np.polyval(np.polyfit(seg.segment, N2, deg=2), seg.segment),
)
ξ = (N2 - seg.N2fit) / seg.N2fit.mean()
N = np.sqrt(seg.N2fit.mean())

_, _, Ptot, m = mixsea.helpers.psd(ξ.data, dp, ffttype="t", window="hamming")

npts = len(ξ.data)
window = signal.get_window("hann", npts)
window /= np.sqrt(np.sum(window**2) / len(window))
kz, psd = signal.periodogram(
    ξ.data, window="hamming", fs=2 * np.pi / dp, detrend="linear", scaling="density"
)
# psd /= np.sinc(kz * dp / 2 / np.pi) ** 2

strain = mixsea.shearstrain.strain_polynomial_fits(
    -1 * gsw.z_from_p(seg.pres_mean.data - 127 + seg.segment.data, seg.latitude.data),
    seg.ctd_temperature.data,
    seg.ctd_salinity.data,
    seg.longitude.data,
    seg.latitude.data,
    # bin_width=64,
    depth_bin=[seg.pres_mean.data],
    window_size=350,
)


plt.loglog(m, Ptot, "x-")
ax = plt.gca()
imix, _ = mixsea.shearstrain.find_cutoff_wavenumber(Ptot, m, 0.22)
dcpy.plots.linex(m[imix[-1]])
print(
    np.trapz(np.real(Ptot[imix]), m[imix]),
    mixsea.shearstrain.gm_strain_variance(m, imix, np.sqrt(strain[0]["N2mean"]))[0],
)

plt.loglog(kz, psd, ".-")
imix, _ = mixsea.shearstrain.find_cutoff_wavenumber(psd, kz, 0.22)
dcpy.plots.linex(kz[imix[-1]], color="C1")
print(
    np.trapz(psd[imix], kz[imix]),
    mixsea.shearstrain.gm_strain_variance(kz, imix, N.data)[0],
)
plt.ylim([1e-3, None])


plt.figure()
plt.plot(strain[0]["N2"])
plt.plot(N2)
plt.plot(strain[0]["N2smooth"], color="C0")
plt.plot(seg.N2fit, color="C1")

plt.figure()
plt.plot(strain[0]["strain"])
plt.plot(ξ)

0.19514013600229327 0.11266552715618153
0.1987434107341453 0.0940107853774018

[<matplotlib.lines.Line2D at 0x7f0b547b8af0>]

_images/1bfa925cfe89f2b63919c6138944d135d294e6e14073aee38919c013a6827f91.png

_images/725718c43ba7a564a0b9887b0c6b8388de771bdf8c89236504ff42520a9fb9c2.png

_images/891f542a4c4c56613039d241a0d65506ffa550a3bf6e201414220c1ca35c7a33.png

GM spectrum is right

# GM 75; Kunze et al (2017)
E0 = 6.3e-5  # nondim spectral energy level
jstar = 3  # peak mode numebr
b = 1300  # m; stratification length scale
N0 = 5.24e-3  # rad/s
K_0 = 5e-6  # m²/s for a mixing efficiency of Γ=0.2
N = N0

kz = 2 * np.pi * np.logspace(-8, -1, 64)  # np.fft.fftfreq(n=1024, d=0.1)[:64]
kzstar = (np.pi * jstar / b) * (N / N0)
ξgm = np.pi * E0 * b / 2 * jstar * kz**2 / (kz + kzstar) ** 2
plt.loglog(kz, ξgm)
print(np.trapz(ξgm, kz))

Sgm, Pgm = mixsea.shearstrain.gm_strain_variance(kz, np.arange(len(kz)), 5.24e-3)
print(Sgm)
plt.plot(kz, Pgm, ".")

0.22007432494266535
0.22007432494266535

[<matplotlib.lines.Line2D at 0x7f0b66d978e0>]

_images/234d34933a77f26afa3249a7eff1e75159d0a0528af5424cd4a97af7d852029e.png

coeffs = segmented.N2.polyfit(deg=2, dim="segment")
smoothed = xr.polyval(segmented.segment, coeffs).rename(
    {"polyfit_coefficients": "N2smooth"}
)
stacked = (
    smoothed.transpose(..., "segment")
    .assign_coords(pres=segmented.pres)
    .stack({"p": ("pressure", "segment")})
)
p0, p1 = stacked.pressure.data[0], stacked.pressure.data[-1]
N2smth = (
    stacked.drop("p").rename({"pres": "pressure"}).swap_dims({"p": "pressure"}).N2smooth
)
N2smth.plot(y="pressure")
mod.N2smooth.plot(y="depth", hue="kind")

[<matplotlib.lines.Line2D at 0x7f0b66b6f8b0>,
 <matplotlib.lines.Line2D at 0x7f0b66b6f8e0>]

_images/1e3fa9fcfeaf582f5ec5e2238688d32924832d1c9d10dc627de11f3dc5da239c.png

Finescale Mixing Parametrizations: Development

Contents

Finescale Mixing Parametrizations: Development#

Check against Kunze & mixsea#

Compare N² estimates:#

Compare to mixsea & Sensitivity studies#