miχPODS: New baseline#

TODO#

  • Add diagnostic for convective “adjustment” or diff == 1

  • fix EUC coordinate shear

  • fix daily composites.

  • finish the deepening metric

  • Add EQUIX, LES

  • write out data reading

Figures#

  • Figure 1 : Jochum background mixing map

  • Figure 2: spectra?

  • mean z-space: u, T, S2?N2?

  • mean EUC space: ν, K, u, T, Ri

  • daily : deepening

  • some daily composite: vertical profile or timeseries

  • seasonal mean heat flux profile.

  • seasonal: monthly mean Jq, SST timeseries

  • ENSO: N2-S2 ENSO + ε-Ri

  • appendix: RigT distribution

Setup#

Hide code cell source
%load_ext autoreload
%load_ext rich
%load_ext watermark

%xmode minimal

import os
import warnings

import cf_xarray as cfxr
import datatree
import dcpy.datatree  # noqa
import distributed
import holoviews as hv
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tqdm
from datatree import DataTree

import xarray as xr

%aimport pump
from holoviews.plotting.bokeh.util import select_legends

from pump import mixpods  # noqa

hv.notebook_extension("bokeh")

cfxr.set_options(warn_on_missing_variables=False)
xr.set_options(keep_attrs=True, display_expand_data=False)

plt.style.use("bmh")
plt.rcParams["figure.dpi"] = 140


%watermark -iv
Hide code cell output
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
Exception reporting mode: Minimal
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.
  from distributed.utils import tmpfile
json       : 2.0.9
matplotlib : 3.7.1
numpy      : 1.24.4
pandas     : 2.0.3
cf_xarray  : 0.8.4
holoviews  : 1.17.0
xarray     : 2023.7.0
pump       : 1.0+273.g892024e.dirty
datatree   : 0.0.12
dcpy       : 0.1.dev387+gd06c937
distributed: 2023.8.0
sys        : 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:40:32) [GCC 12.3.0]
tqdm       : 4.66.1
Hide code cell source
if "client" in locals():
    client.close()  # noqa
    del client  # noqa
if "cluster" in locals():
    cluster.close()  # noqa

import ncar_jobqueue

cluster = ncar_jobqueue.NCARCluster(
    local_directory="/local_scratch/pbs.$PBS_JOBID/dask/spill",
)
cluster.adapt(minimum_jobs=2, maximum_jobs=8)
client = distributed.Client(cluster)
client
Hide code cell output

Client

Client-03546858-3d3b-11ee-a985-3cecef1ac722

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/dcherian/casper/proxy/8787/status

Cluster Info

Introduction#

Analysis of χpod shear mixing measurements recorded over approximately 15 years has revealed that mixing varies on multiple timescales, from hourly to interannual, in the tropical Pacific and Indian Oceans (Moum, 2021; Cherian et al., 2020b; Warner et al., 2016; Warner and Moum, 2019; Pujiana et al., 2018; Moum et al., 2009, 2013). These results show that equatorial mixing is organized in large-scale coherent patterns.

We pask whether these patterns, a decade-long observed mixing dataset (from χpods) to develop mixing process-oriented diagnostics (miχPODs) that assess the variability of parameterized equatorial mixing from monthly to decadal timescales; and its accumu- lated impact on model biases by applying miχPODs to a suite of existing CMIP6-class forced ocean and coupled model experiments from NCAR and GFDL.

Following Large and Gent (1999), the comparison of a 1D mixing model with LES or DNS is clean, because the forcing is the same, and direct, because the evaluation compares turbulence quantities (despite shortcomings such as model errors and idealized forcing). In contrast, the performance of a mixing scheme in global and regional model simulations is tested by comparing the simulated mean state to an observed mean state after a long integration (for e.g. Blanke and Delecluse, 1993; Robitaille and Weaver, 1995; Large and Gent, 1999; Gutjahr et al., 2020). This comparison is not clean, because of errors in forcing fields, nor is it direct, since mean state properties are not a direct output of the mixing scheme.

We propose diagnostics that evaluate the representation of the variability in parameterized equatorial mixing using the multiyear χpod turbulence dataset. While this comparison is not clean, it is direct, a substantial improvement to the current state-of-the-art approach.

These diagnostics are direct evaluations of the models’ ability to simulate the observed variability of ocean vertical mixing and reproduce large-scale patterns recorded in equatorial mixing observations.

Datasets#

Observations#

  1. TAO \(T, S, u, v\) and turbulence (Moum et al XXX)

  • Note which depth χpods are being used

  • Add a figure showing sampling and gaps

Model Simulations#

The effect of these parameter changes are examined in a IPCC-class climate model: the development version of Community Earth System Model version 3 (CESM3), using the Modular Ocean Model version 6 (MOM6) ocean component, and the K Profile Paramerization vertical mixing scheme (Large et al., 1994).

Ocean-ice simulations (“g case”) are integrated using JRA-55 forcing dataset.

Viscosity and diffusivity in the model#

  • TODO write down KPP and Jochum equations

Total diffusivity \(K_D^{Total} = K_D + K_D^{Jochum}(x, y)+ K_D^{KPP}(x,y,z,t) \) . The constant \(K_D\) and the spatially varying \(K_D^{Jochum}(x,y)\) model the diffusivity due to breaking internal waves. \(K_D^{KPP}\) models surface layer mixing and mixing due to shear instability.

  • Split \(K_D^{KPP}\) too.

The same decomposition applies to viscosity \(K_V\).

staticfile = "/glade/campaign/cgd/oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb/run/*static*.nc"
static = xr.open_mfdataset(staticfile).squeeze()
(
    static.Kd_bkgnd.cf.sel(Z=200, method="nearest").hvplot.quadmesh(
        cmap="spectral_r",
        cnorm="log",
        clim=(5e-7, 2e-4),
        title="Jochum (2009) background diffusivity",
    )
    + static.Kv_bkgnd.cf.sel(Z=200, method="nearest").hvplot.quadmesh(
        cmap="spectral_r",
        cnorm="log",
        clim=(5e-7, 2e-4),
        title="Jochum (2009) background viscosity",
    )
).cols(1)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  layout_plot = gridplot(

Parameters#

We perturb the following parameters

  1. “Background mixing” parameters used to simulate mixing by unresolved processes such as breaking internal waves.

    1. A spatially varying background diffusivity \(K_D^{Jochum}(x,y)\) using the formulation by Jochum (2003)

    2. A spatially varying background viscosity \(K_V^{Jochum}\) derived from \(K_D^{Jochum}\) using a Prandtl number of 5 (CHECK).

    3. A constant background viscosity \(K_V\) added to (B).

  2. KPP Boundary layer scheme:

    1. Bulk Richardson number criterion \(Ri_c\) which controls the depth of the actively mixing layer.

  3. KPP shear mixing scheme:

    1. Critical Richardson number \(Ri_0\) that controls when mixing occurs.

    2. Maximum shear mixing viscosity \(ν_0\)

Presented Configurations#

Choices for the parameter changes were influenced by experiments matching 1D KPP simulations to the Large Eddy Simulations of Whitt et al (2022). The range of parameter values explored here is not exhaustive and is only intended to demonstrate the utility of these metrics in tuning and assessing model performance.

We run baseline cases for a full cycle 1958-2018 (inclusive).

  1. “old” baseline with KD=1e-5, KV=2e-4 (NCAR/MOM6#238)

  2. new baseline : KD=0, KV=0

Parameter sensitivities are explored by branching the new simulation starting at 2003/01/01 using the baseline simulation. We anticipate quick adjustments to the new parameters.

  1. old baseline KD=1e-5, KV=2e-4, kpp.lmd.004 with KPP ν0=2.5e-3, Ric=0.2, Ri0=0.5

  2. new baseline KD=0, KV=0, kpp.lmd.004 with KPP ν0=2.5e-3, Ric=0.2, Ri0=0.5

  3. new baseline KD=0, KV=0, kpp.lmd.005 with KPP ν0=2.5e-3, Ric=0.3, Ri0=0.5

For all simulations, we save a large number of variables at the TAO mooring locations every time step, which is an hour.

This paper primarily describes output at (0, 140W).

Define quantities#

Shear is estimated using central differences.

Summary#

  1. Turning down the background visc by a factor of 40, (2e-4 → 5e-5 m2/s)

    1. sharpens the EUC relative to TAO

    2. Some background visc is necesary below the EUC to match known physics (Peters et al, 2005)

  2. Using Ri_c=0.2, so shallower KPP surface layer, is key as Dan mentioned.

  3. Median Ri / “marginal stability’ is not very useful diagnostic at 140W

Read#

MOM6#

Look at catalog#

The catalog only contains surface variables, full depth variables, NO SECTIONS YET.

I use it to keep track of available cases

catalog = pump.catalog.open_mom6_catalog()
catalog.df["shortname"].unique()

array(['baseline', 'kpp.lmd.004', 'baseline.001', 'baseline.epbl.001',
       'baseline.hb', 'baseline.kpp.lmd.002', 'baseline.kpp.lmd.003',
       'baseline.kpp.lmd.004', 'new_baseline.hb',
       'new_baseline.kpp.lmd.004', 'new_baseline.kpp.lmd.005'],
      dtype=object)
catalog.df

casename stream path baseline levels frequency variables shortname description
0 gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 150 N/A [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline baseline N=150
1 gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 150 daily [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline baseline N=150
2 gmom.e23.GJRAv3.TL319_t061_zstar_N150.kpp.lmd.... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 150 N/A [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5, N=150
3 gmom.e23.GJRAv3.TL319_t061_zstar_N150.kpp.lmd.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 150 daily [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5, N=150
4 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 N/A [Heat_PmE, KE, N2_int, Rd_dx, SSH, SSU, SSV, T... baseline.001 baseline
5 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... baseline.001 baseline
6 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [Heat_PmE, KE, Rd_dx, SSH, SSU, SSV, T_adx_2d,... baseline.001 baseline
7 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [Rd_dx, SSH, SSU, SSV, mass_wt, mlotst, oml, o... baseline.001 baseline
8 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 N/A [SSH, SSU, SSV, T_advection_xy, T_lbdxy_cont_t... baseline.epbl.001 ePBL
9 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.epbl.001 ePBL
10 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, sos, speed, tos] baseline.epbl.001 ePBL
11 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [T_advection_xy, T_lbdxy_cont_tendency, Th_ten... baseline.epbl.001 ePBL
12 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 N/A [Heat_PmE, KE, KPP_NLT_temp_budget, N2_int, Rd... baseline.hb baseline
13 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... baseline.hb baseline
14 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [Heat_PmE, KE, Rd_dx, SSH, SSU, SSV, SW_pen, T... baseline.hb baseline
15 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 daily [Rd_dx, SSH, SSU, SSV, mass_wt, mlotst, oml, o... baseline.hb baseline
16 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.hb baseline
17 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 N/A [KPP_NLT_temp_budget, SSH, SSU, SSV, T_advecti... baseline.kpp.lmd.002 KPP Ri0=0.5
18 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.kpp.lmd.002 KPP Ri0=0.5
19 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline.kpp.lmd.002 KPP Ri0=0.5
20 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.kpp.lmd.002 KPP Ri0=0.5
21 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 N/A [KPP_NLT_temp_budget, SSH, SSU, SSV, T_advecti... baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
22 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
23 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
24 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.kpp.lmd.003 KPP Ri0=0.5, Ric=0.2,
25 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 N/A [KPP_NLT_temp_budget, SSH, SSU, SSV, T_advecti... baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
26 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [h, so, thetao, uhGM, uhml, umo, uo, vhGM, vhm... baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
27 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [SSH, SSU, SSV, mlotst, oml, sos, speed, tos] baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
28 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.... wci /glade/campaign/cgd/oce/projects/pump/cesm/gmo... old 65 monthly [KPP_NLT_temp_budget, T_advection_xy, T_lbdxy_... baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
29 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 N/A [Heat_PmE, KE, KPP_NLT_temp_budget, N2_int, Rd... new_baseline.hb KD=0, KV=0
30 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... new_baseline.hb KD=0, KV=0
31 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [Heat_PmE, KE, KPP_NLT_temp_budget, Rd_dx, SSH... new_baseline.hb KD=0, KV=0
32 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 daily [Rd_dx, SSH, SSU, SSV, mass_wt, mlotst, oml, o... new_baseline.hb KD=0, KV=0
33 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 N/A [Heat_PmE, KE, KPP_NLT_temp_budget, N2_int, Rd... new_baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
34 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... new_baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
35 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [Heat_PmE, KE, KPP_NLT_temp_budget, Rd_dx, SSH... new_baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
36 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 daily [Rd_dx, SSH, SSU, SSV, mass_wt, mlotst, oml, o... new_baseline.kpp.lmd.004 KPP ν0=2.5, Ric=0.2, Ri0=0.5
37 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... combined /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 N/A [Heat_PmE, KE, KPP_NLT_temp_budget, N2_int, Rd... new_baseline.kpp.lmd.005 KPP ν0=2.5, Ri0=0.5
38 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... h /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [N2_int, agessc, h, rhopot0, so, thetao, uhGM,... new_baseline.kpp.lmd.005 KPP ν0=2.5, Ri0=0.5
39 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... hm /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 monthly [Heat_PmE, KE, KPP_NLT_temp_budget, Rd_dx, SSH... new_baseline.kpp.lmd.005 KPP ν0=2.5, Ri0=0.5
40 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_basel... sfc /glade/campaign/cgd/oce/projects/pump/cesm/gmo... new 65 daily [Rd_dx, SSH, SSU, SSV, mass_wt, mlotst, oml, o... new_baseline.kpp.lmd.005 KPP ν0=2.5, Ri0=0.5

Pick these simulations

SIMS = [
    "baseline.hb",
    "baseline.kpp.lmd.004",
    "new_baseline.hb",
    "new_baseline.kpp.lmd.004",
    "new_baseline.kpp.lmd.005",
]
catalog_sims = catalog.search(shortname=SIMS)

Subset catalog#

Note

Use baseline.hb instead of baseline.001. This was run by Anna-Lena, and does not have gaps, and has heat budget output.

catalog_sub = catalog_sims.search(stream="combined")
with pd.option_context("display.max_colwidth", None):
    display(catalog_sub.df[["shortname", "casename", "path"]])  # noqa

shortname casename path
0 baseline.hb gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb /glade/campaign/cgd/oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.hb/run/jsons/combined.json
1 baseline.kpp.lmd.004 gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods /glade/campaign/cgd/oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods/run/jsons/combined.json
2 new_baseline.hb gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb /glade/campaign/cgd/oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb/run/jsons/combined.json
3 new_baseline.kpp.lmd.004 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods /glade/campaign/cgd/oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods/run/jsons/combined.json
4 new_baseline.kpp.lmd.005 gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods /glade/campaign/cgd/oce/projects/pump/cesm/gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods/run/jsons/combined.json
backup#

I tried catalog.search but that is ugly and didn’t help with reading SST etc. later

catalog_sub = pd.concat(
    [
        cat.df
        for cat in [
            # pick two "old" baseline simulations
            catalog.search(
                stream="combined",
                levels=65,
                shortname=["baseline.hb", "baseline.kpp.lmd.004"],
            ),
            # and all new baseline simulations
            catalog.search(stream="combined", levels=65, baseline="new"),
        ]
    ]
)

with pd.option_context("display.max_colwidth", None):
    display(catalog_sub[["shortname", "casename", "path"]])  # noqa

This is what I was doing before switching to catalog.search

catalog_sub = {
    k: catalog[k]
    for k in catalog.keys()
    if k == "kpp.lmd.004" or ("baseline" in k and "150" not in k)
}
display(catalog_sub)
{
    'baseline': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods'),
    'kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.hb': ('KD=0, KV=0', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb'),
    'new_baseline.kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.kpp.lmd.005': (
        'KPP ν0=2.5, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods'
    )
}

load subset to tree#

Since there are no sections in the catalog, and we do a bunch of custom things, loop over entries and load the section output.

%autoreload

datasets = {}
with warnings.catch_warnings():
    warnings.simplefilter("ignore", category=UserWarning)
    warnings.simplefilter("ignore", category=FutureWarning)
    for _, (casename, shortname, description) in tqdm.tqdm(
        catalog_sub.df[["casename", "shortname", "description"]].iterrows()
    ):
        datasets.update(
            {
                shortname: mixpods.load_mom6_sections(casename).assign_attrs(
                    {"title": description}
                )
            }
        )
5it [00:10,  2.12s/it]
tree = DataTree.from_dict(datasets)
tree

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

Add TAO#

tao_gridded = mixpods.load_tao()
tree["TAO"] = DataTree(tao_gridded)
tree = tree.dc.reorder_nodes(["TAO", ...])
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:265: UserWarning: The specified chunks separate the stored chunks along dimension "depth" starting at index 42. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:265: UserWarning: The specified chunks separate the stored chunks along dimension "time" starting at index 199726. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:265: UserWarning: The specified chunks separate the stored chunks along dimension "longitude" starting at index 2. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(

Add LES, microstructure#

Check SST time series#

sst = (
    xr.Dataset(
        {
            casename: node["sst"].reset_coords(drop=True)
            for casename, node in tree.children.items()
        }
    )
    .to_array("case")
    .load()
)
sst.rename("sst").resample(time="D").mean().hvplot.line(muted_alpha=0, by="case")

Post-process catalog subset#

Slice#

tree = tree.sel(time=slice("2003", "2018"))
tree

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

N2-S2 histograms#

ref = tree["TAO"].ds.reset_coords(drop=True).cf.sel(Z=slice(-120, None))
counts = np.minimum(ref["S2"].cf.count("Z"), ref["N2T"].cf.count("Z")).load()


def calc_histograms(ds):
    ds = ds.copy()
    ds["tao_mask"] = counts.reindex(time=ds.time, method="nearest") > 5
    ds["tao_mask"].attrs = {
        "description": "True when there are more than 5 5-m T, u, v in TAO dataset"
    }
    # I did try masking out the model like the data
    # ds = ds.where(ds.tao_mask)
    return ds.update(mixpods.pdf_N2S2(ds))


tree = tree.map_over_subtree(calc_histograms)

daily composite#

dailies = tree.map_over_subtree(mixpods.daily_composites)
dailies
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

euc relative frame#

Calculate#

This one is a bit intensive

cluster.adapt(minimum_jobs=10)

<distributed.deploy.adaptive.Adaptive object at 0x2afadd620160>
# could be cleaned up
euc = mixpods.bin_to_euc_centered_coordinate(tree)
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
# merge tree
euc_mean = mixpods.average_euc(euc)

euc = euc.dc.insert_as_subtree("mean", euc_mean)
euc.to_zarr("mom6-euc-coordinate.zarr")
cluster.adapt(minimum_jobs=2)

<distributed.deploy.adaptive.Adaptive object at 0x2afb355fb1c0>
Read existing cache#
euc = datatree.open_datatree("mom6-euc-coordinate.zarr", engine="zarr")
euc_mean = euc.dc.extract_leaf("mean")
# tree = tree.dc.insert_as_subtree("euc", euc)
Merge in microstructure#
micro_zeuc = datatree.open_datatree(
    os.path.expanduser("~/datasets/microstructure/equix-tiwe-zeuc.average.nc")
).load()
for nodename, node in micro_zeuc.children.items():
    node["u"].attrs["standard_name"] = "sea_water_x_velocity"
    node["v"].attrs["standard_name"] = "sea_water_y_velocity"
    node["KT"].attrs["standard_name"] = "ocean_vertical_heat_diffusivity"
    node["ν"].attrs["standard_name"] = "ocean_vertical_momentum_diffusivity"

    del node["Rig_T"].attrs["standard_name"]
    del node["N2T"].attrs["standard_name"]
    # euc_mean[nodename] = node

for nodename, node in micro_zeuc.children.items():
    euc_mean[nodename] = node

Validate before continuing#

mixpods.validate_tree(tree)

Persist tree#

tree = mixpods.persist_tree(tree)
euc_mean.load()

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

Scales of comparison#

Comparing a coarse climate model at 2/3° lateral grid spacing, and a time step of an hour, to microstructure observations that capture temperature variability at 7-15Hz is challenging.

The equator is a special place with strong diurnal variability in turbulence.

Capturing the daily timescale at the equator has been a vital benchmark for the KPP scheme (Large and Gent, 1999).

Interestingly we find that all simulations here reproduce a qualitatively similar variability at the daily frequency, even though variance is depressed at sub-daily and super-daily frequencies.

Conclusion:

  1. We should be comparing at daily frequency, not hourly.

  2. Vertical wavenumber spectra are equivalent to TAO

TODO: Frequency spectra comparison#

  • Now obvious difference between clockwise and counter-clockwise so I’ve chosen to just do the spectrum of the magnitude

  • calculate shear before regridding

  • There is sensitivity to coordinate location chosen.

  • I picked one year : 2010

  • Models are doing quite well at daily frequency at the 40-50m range.

f, ax = plt.subplots(1, 2, sharey=True)
for axis, t, depth in zip(ax, [tree, euc], [-45, 30]):
    for name, node in t.children.items():
        da = node.ds["wz"].reset_coords(drop=True)
        if "zeuc" in da.dims:
            da.zeuc.attrs["axis"] = "Z"
        dcpy.ts.PlotSpectrum(
            np.abs(da.cf.sel(Z=depth, method="nearest").sel(time="2010")),
            multitaper=False,
            nsmooth=12,
            label=name,
            lw=0.75,
            ax=axis,
        )
ax[0].set_ylim([1e-9, None])
ax[0].legend()
ax[0].set_title("Shear spectra, depth=-40m")
ax[1].set_title("Shear spectra, zeuc=+30m")
f.set_size_inches((10, 4))

_images/9608b031c56f2de5e077a6e7bb38f0a102a0894223a764b3ac1d49b0b39bb635.png

Vertical wavenumber spectra#

  • Add EQUIX, LES

TAO drops at 0.04 m^-1 compared to EQUIX, LES.

new_baseline.kpp.lmd.00X drop off at the same point.

Consider tweaked model and TAO to be equivalent?

Choosing to do w.cf.differentiate("Z") . This will use central differences on all dataset. wz uses centered differences for TAO and forward differences for models. Yields different results.

%autoreload

spec_z = mixpods.hvplot_spectra(
    tree.map_over_subtree(
        lambda node: (
            node["w"]
            .cf.differentiate("Z")
            .cf.sel(Z=slice(-200, 0))
            .cf.interp(Z=np.arange(-105, -52, 5))
        )
    ),
    varname="w",
    dim="Z",
).opts(title="Z=slice(-100, -50)")
spec_z
Task exception was never retrieved
future: <Task finished name='Task-711064' coro=<Client._gather.<locals>.wait() done, defined at /glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/distributed/client.py:2232> exception=AllExit()>
Traceback (most recent call last):
  File "/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/distributed/client.py", line 2241, in wait
    raise AllExit()
distributed.client.AllExit

Mean#

Mean State Variable Profiles#

  1. A lot more shear, \(S, S^2\) just above the EUC when visc is turned down!

  2. We are very slightly lower on \(S^2\), \(N_T^2\) in the top 80m, compare TAO, kpp.lmd.004, new_baseline.kpp.lmd.004.

  3. Using the standard Ri_c=0.3, so deeper KPP surface layer, decreases \(S^2\), \(N_T^2\) in the top 60m. new_baseline.kpp.lmd.004 vs new_baseline.kpp.lmd.005

for context, Peters et al (1995) is interesting:

  1. Variability patterns at 50-350 m are distinctly different from the upper 50 m containing the diurnal cycle of mixing

  2. Large-scale shear of wavenumbers k < 0.01 cpm is associated with the slowly varying EUC and EIC. Large-scale is the dominant component of total shear above 50 m and in the thermostad, 170-270 m.

  3. Fine-scale shear, 0.01cpm < k < 0.5 cpm,provides the dominant component of total rms shear and exceeds the large-scale component in thick layersaroundthe coresof EUC and EIC, where Fr < 1, at 50-170 m and below 270 m.

  4. Variations of fine-scale shear cause variations in turbulent mixing; the large-scaleshear alone is a poor predictor of mixing.

  5. Lacking a model that links large-scale, fine-scale, and turbulent flow components,our service to equatorial modelers consists of describing general levels of eddy diffusivities and their variability patterns.

These models should not be resolving “finescale” shear, and the mixing is not well correlated with “large-scale shear”, so we need a “background” diffusivity/viscosity to make things work.

But see that std(\(S^2\)) is a lot stronger above the EUC core, relative to TAO for new_baseline.kpp.lmd.004 and new_baseline.kpp.lmd.005

Conclusion Tweaking KPP works nicely in the deep-cycle layer. We see sharper top of the EUC. However bottom half is too sharp.

S2 = mixpods.plot_profile_fill(tree, "S2", "S^2")
N2 = mixpods.plot_profile_fill(tree, "N2T", "N_T^2")
u = mixpods.plot_profile_fill(tree, "sea_water_x_velocity", "u")
T = mixpods.plot_profile_fill(tree, "sea_water_potential_temperature", "T")
# Ri = mixpods.plot_median_Ri(tree)
%autoreload

plot = (
    (S2 + N2 + u + T)  # + Ri
    .cols(5)
    .opts(
        hv.opts.Curve(xlim=(-350, 0), frame_height=600),
        hv.opts.Layout(toolbar="left"),
        hv.opts.Overlay(legend_offset=(0, 200)),
        clone=True,
    )
)
select_legends(plot, figure_index=4, legend_position="right")
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(

Mean profiles in EUC-relative coordinate#

tao_chi_counts = euc["TAO/chi"].count("time").load()

euc_mean["TAO"] = euc_mean["TAO"].update(
    euc_mean["TAO"]
    .ds.cf[
        [
            "chi",
            "eps",
            "ocean_vertical_heat_diffusivity",
            "ocean_vertical_momentum_diffusivity",
        ]
    ]
    .where(tao_chi_counts > 2000)
)
tao_chi_counts.hvplot.line(title="TAO χ number of hourly obs in bin")

h = {
    varname: mixpods.map_hvplot(
        lambda node, name, muted: (
            node.ds.cf[varname]
            .reset_coords(drop=True)
            .hvplot.line(
                ylabel=varname,
                label=name,
                logx=varname != "sea_water_x_velocity",
                invert=True,
            )
        ),
        euc_mean,
    )
    for varname in [
        "sea_water_x_velocity",
        "sea_water_potential_temperature",
        "chi",
        "eps",
        "ocean_vertical_heat_diffusivity",
        "ocean_vertical_momentum_diffusivity",
    ]
}
h["chi"].opts(ylim=(1e-9, 1e-4))
h["eps"].opts(ylim=(1e-9, 1e-4))
h["ocean_vertical_heat_diffusivity"].opts(ylim=(5e-7, 3))
h["ocean_vertical_momentum_diffusivity"].opts(ylim=(1e-6, 1e-1))

h2 = {
    varname: mixpods.map_hvplot(
        lambda node, name, muted: node.ds[varname]
        .reset_coords(drop=True)
        .hvplot.line(ylabel=varname, label=name, logx=False, invert=True),
        euc_mean,
    )
    for varname in ["Rig_T"]
}

h

plot = (
    hv.Layout(list(h.values()) + [h2["Rig_T"].opts(ylim=(0, 3))])
    .opts(
        hv.opts.Curve(frame_width=150, frame_height=300, xlim=(-150, 150)),
        hv.opts.Layout(shared_axes=True),
        hv.opts.Overlay(show_legend=True, show_grid=True, legend_position="right"),
        *mixpods.HV_TOOLS_OPTIONS,
        *mixpods.PRESENTATION_OPTS,
    )
    .cols(4)
)
select_legends(plot, figure_index=3, legend_position="right")
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(

A deeper look at mixing below the EUC#

We find that tuning KPP’s shear scheme to better represent deep cycle turbulence by reducing \(Ri_0\) to 0.5 from 0.8 where median \(Ri_g\) is lower (0.25), prevents it from kicking in in the lower half of the EUC where median \(Ri_g\) is higher (1).

  • viscosity is too low below EUC

  • but diffusivity may be OK?

Microstructure observations from short term experiments (TIWE, EQUIX) do show higher viscosities \(ν \sim\) \SI{1e-4}{\meter\squared\per\second} below the EUC maximum.

We choose to preserve only TAO χpods above 90m depth, so there is very little data below the EUC core. For these data points, bins with less than 2000 hourly measurements are excluded.

With reduced \(Ri_0 \sim 0.5\), the KPP scheme is unable to reproduce these viscosities resulting a sharper EUC.

Peters et al (1995):

  1. Fine-scale shear, 0.01cpm < k < 0.5 cpm,provides the dominant component of total rms shear and exceeds the large-scale component in thick layersaroundthe coresof EUC and EIC, where Fr < 1, at 50-170 m and below 270 m.

  2. Variations of fine-scale shear cause variations in turbulent mixing; the large-scale shear alone is a poor predictor of mixing.

These simulations should not be resolving “finescale” shear, and the mixing is not well correlated with “large-scale shear”, so we need a “background” diffusivity/viscosity to make things work.

Conclusion Indeed we see that additional background viscosity is necessary — compare new_baseline.hb to new_baseline.kpp.lmd.004.

One might consider tweaking shear scheme to do better at deep-cycle, and background viscosity to represent mixing due to finescale shear below the EUC core. This seems physically justified.

Mean synoptic metrics#

Read#

uo = mixpods.read_stream_vars(
    catalog_sims,
    stream="hm",
    variables="uo",
).mean("time")
uo
--> The keys in the returned dictionary of datasets are constructed as follows:
	'shortname.stream'


100.00% [4/4 00:01<00:00]






<xarray.DataArray (dataset: 4, zl: 65, yh: 458, xq: 540)>
dask.array<chunksize=(1, 65, 458, 540), meta=np.ndarray>
Coordinates:
  * xq       (xq) float64 -286.3 -285.7 -285.0 -284.3 ... 71.0 71.67 72.33 73.0
  * yh       (yh) float64 -79.2 -79.08 -78.95 -78.82 ... 87.55 87.64 87.71 87.74
  * zl       (zl) float64 1.25 3.75 6.25 8.75 ... 5.378e+03 5.627e+03 5.876e+03
  * dataset  (dataset) object 'baseline.hb' ... 'new_baseline.kpp.lmd.005'
thetao = mixpods.read_stream_vars(
    catalog_sims,
    stream="hm",
    variables="thetao",
).mean("time")
thetao
--> The keys in the returned dictionary of datasets are constructed as follows:
	'shortname.stream'


100.00% [4/4 00:00<00:00]






<xarray.DataArray (dataset: 4, zl: 65, yh: 458, xh: 540)>
dask.array<chunksize=(1, 65, 458, 540), meta=np.ndarray>
Coordinates:
  * xh       (xh) float64 -286.7 -286.0 -285.3 -284.7 ... 70.67 71.33 72.0 72.67
  * yh       (yh) float64 -79.2 -79.08 -78.95 -78.82 ... 87.55 87.64 87.71 87.74
  * zl       (zl) float64 1.25 3.75 6.25 8.75 ... 5.378e+03 5.627e+03 5.876e+03
  * dataset  (dataset) object 'baseline.hb' ... 'new_baseline.kpp.lmd.005'
sst = mixpods.read_stream_vars(catalog_sims, stream="sfc", variables="tos").mean("time")
sst
--> The keys in the returned dictionary of datasets are constructed as follows:
	'shortname.stream'


100.00% [5/5 00:03<00:00]







<xarray.DataArray (dataset: 5, yh: 458, xh: 540)>
dask.array<chunksize=(1, 458, 540), meta=np.ndarray>
Coordinates:
  * xh       (xh) float64 -286.7 -286.0 -285.3 -284.7 ... 70.67 71.33 72.0 72.67
  * yh       (yh) float64 -79.2 -79.08 -78.95 -78.82 ... 87.55 87.64 87.71 87.74
  * dataset  (dataset) object 'baseline.hb' ... 'new_baseline.kpp.lmd.005'

Mean cold tongue SST contours#

sst.hvplot.contour(
    levels=(23, 24, 25, 25.5, 26, 26.5),
    by="dataset",
    xlim=(-170, -80),
    ylim=(-10, 10),
)

sst.hvplot.contour(
    levels=(23, 24, 25, 25.5, 26, 26.5),
    groupby="dataset",
    xlim=(-170, -80),
    ylim=(-10, 10),
)

Along-equator section#

  • Lines of depth of thermocline and depth of EUC along equator

Cross-equator section :#

  1. integral(u dz),

  2. depth of dT/dz max at each latitude

The daily timescale#

Type in motivation

Conclusions ????

Daily cycle without wind dependencies#

  • a profile?

Daily composites: ε (Moum et al, 2023)#

This figure is hard to interpret!

I think a lot of this is bias in KPP surface layer vs actively mixing layer in obs. We can write better diagnostics to check this (Moum et al, 2023)

  1. The 89m χpod comparison is quite interesting. Suggests we do need more background visc

plot_kt = mixpods.plot_daily_composites(
    dailies, ["ocean_vertical_heat_diffusivity"], logy=True
)
plot_kt
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:588: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  plot = gridplot(plots[::-1],
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:588: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  plot = gridplot(plots[::-1],
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:647: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  plot = gridplot([r1, r2])
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:647: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  plot = gridplot([r1, r2])

%autoreload

plot_eps = mixpods.plot_daily_composites(dailies, ["eps"], logy=True)
# plot_kt = mixpods.plot_daily_composites(dailies, ["ocean_vertical_heat_diffusivity"], logy=True)
(plot_eps).opts(toolbar="left")
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:588: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  plot = gridplot(plots[::-1],
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:588: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  plot = gridplot(plots[::-1],
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:647: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  plot = gridplot([r1, r2])
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:647: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  plot = gridplot([r1, r2])

%autoreload

mixpods.plot_daily_composites(
    dailies, ["S2", "N2", "Rig_T"], logy=False, legend=False
).opts(hv.opts.GridSpace(show_legend=True, width=200), hv.opts.Overlay(frame_width=200))

Deepening of the boundary layer#

  • TODO: map to local time, use shortwave as reference

We follow the Moum et al (2023) definition

We quantify the progression in time of the deepening of DC turbulence shown in Figures 7–9 by finding the time of the maximum rate of increase of \(𝜖\), \(𝜖_{𝑚𝑎𝑥}^𝑡\) , from each series after filtering each series at 1 hour. This analysis yields the deepening trajectories shown in Figures 10a-c, at the depths where an unambiguous determination of \(𝜖^{𝑚𝑎𝑥}_𝑡\) could be made. These trajectories show the arrival of DC turbulence to be progressively delayed at all depths with weakening winds.

eps = dailies.dc.subset_nodes("eps").dc.concatenate_nodes().eps
idx = (
    eps.diff("hour")  # .isel(depth=0)
    .rolling(hour=5, center=True, min_periods=1)
    .mean()
    .idxmax("hour")
    .assign_attrs(long_name="time of max $ε_t$")
)
# mask = ((-1 * idx.diff("depth", label="lower")) < 0).reindex(
#    depth=idx.depth, fill_value=False
# )
# idx = xr.where(mask, idx + 12, idx)
start = (3, -30)
end = 12

plot = hv.Layout(
    [
        da.hvplot.line(y="depth", by="node")
        * hv.Curve([start, (end, -30 - 6 * (end - start[0]))]).opts(
            line_dash="dashed", line_color="k"
        )
        for _, da in idx.groupby("tau_bins")
    ]
).opts(
    hv.opts.Curve(frame_height=300, frame_width=250),
    *mixpods.HV_TOOLS_OPTIONS,
    *mixpods.PRESENTATION_OPTS,
)
select_legends(plot, 2, legend_position="right")
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(

Seasonal timescale (Moum et al, 2013)#

Conclusion :

  • Timeseries : Jq is not too bad. Potentially some phase offset at April/May. Not sure Why

  • Profiles : Jq biases are significant.

Seasonal mean profiles#

vertical = tree.map_over_subtree(
    lambda node: (
        node.cf[["Jq", "ocean_vertical_heat_diffusivity"]]
        .cf.rename({"ocean_vertical_heat_diffusivity": "KT"})
        .cf.sel(Z=slice(-120, -20))
        .groupby("time.season")
        .mean()
        # reorder
        .sel(season=["DJF", "MAM", "JJA", "SON"])
        .load()
    )
)

seasonal = {}
seasonal["Jq"] = mixpods.map_hvplot(
    lambda ds, name, muted: (
        ds["Jq"]
        .cf.rename({"Z": "depth"})
        .hvplot.line(y="depth", col="season", label=name, muted=muted)
    ),
    vertical.children,
)
seasonal["KT"] = mixpods.map_hvplot(
    lambda ds, name, muted: (
        ds["KT"]
        .cf.rename({"Z": "depth"})
        .hvplot.line(logx=True, y="depth", col="season", label=name, muted=muted)
    ),
    vertical.children,
)
proc = vertical.map_over_subtree(lambda node: node.cf.rename({"Z": "Z"}))
layout = mixpods.hvplot_facet(
    proc,
    varname="Jq",
    col="season",
)
layout.opts(
    hv.opts.Curve(invert_axes=True, xrotation=20, frame_width=160, frame_height=300),
    hv.opts.Overlay(shared_axes=True, labelled=["x"]),
    *mixpods.HV_TOOLS_OPTIONS,
    *mixpods.PRESENTATION_OPTS,
)
select_legends(layout, figure_index=3, legend_position="right")
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(

Monthly mean timeseries Jq, SST#

  • Note this is sensitive to order over averaging. I think it’s sane to create monthly climatology for each depth and then average across depth.

  • Not sure TAO “SST” here is the best. Let’s use another product.

sub = tree.map_over_subtree(
    lambda node: (
        node.cf[["Jq", "sea_surface_temperature", "ocean_vertical_heat_diffusivity"]]
        .pipe(np.abs)
        .cf.rename(
            {"sea_surface_temperature": "sst", "ocean_vertical_heat_diffusivity": "KT"}
        )
        .cf.sel(Z=slice(-60, -20))
        .groupby("time.month")
        .mean()
        .cf.mean("Z")
        .load()
    )
).dc.concatenate_nodes()
h = (
    (sub.Jq.hvplot.line(by="node") + sub.sst.hvplot.line(by="node", legend=False))
    .opts(
        hv.opts.Curve(frame_width=500),
        *mixpods.PRESENTATION_OPTS,
        *mixpods.HV_TOOLS_OPTIONS,
    )
    .cols(1)
)
h
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  layout_plot = gridplot(

ENSO Timescale#

Conclusion: Models do quite well! Tweaking shows us better at marginal stability and reproducing higher magnitude S2, N2. Dynamic range between ENSO phases seems OK. Some reproduction of variability in ENSO phase in ε-Ri plot.

Stability Diagram#

  1. ~Major change in La-Nina Warming for the new-baseline runs. We need to check ONI closely.~

    • fixed by just using the obs ONI everywhere.

      • This is probably OK since it shouldn’t drift too far away.

      • And we now average over the same time periods.

    • It is a pain to estimate this with branch runs,

    • but could be done if we spliced in the SST from the baseline case

  2. I still can’t reproduce the Warner and Moum (2019) figure where the El-Nino warming phase overlaps a lot less with La-Nina cooling

  • I’m only using data in the “deep cycle layer”.

    • have not matched vertical grid spacing yet.

    • could interpolate to TAO χpod depths

  • think about resampling to daily. Currently hourly.

  • think about instrumental error on vertical gradients for S2, N2;

    • if added to model, would it change the bottom tail.

%autoreload

mixpods.plot_stability_diagram_by_dataset(tree, nrows=2)

_images/9906b7b0d2696f5da29fa2faf152ee3edfa3076d27400ddce1d1a1642056e405.png

ε-Ri histograms#

%autoreload

mixpods.map_hvplot(
    lambda dt, name, muted: mixpods.plot_eps_ri_hist(
        dt["eps_ri"].load(), label=name, muted=muted
    ),
    tree.children,
).opts(
    hv.opts.Curve(ylim=(1e-8, 3e-6)),
    hv.opts.GridSpace(show_legend=True),
    hv.opts.Overlay(show_grid=True, legend_position="bottom_left"),
    *mixpods.HV_TOOLS_OPTIONS,
)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:588: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  plot = gridplot(plots[::-1],

Discussion#

Process based metrics point a way forward to tuning of a model vertical mixing scheme.

We also find structural deficiencies still exist.

Supplementary Material#

Turbulence Variable distributions#

  • Again hourly vs daily

Maybe just talk about Ri distributions. don’t think the others are useful.

handles = [
    mixpods.plot_distributions(tree, "chi", bins=np.linspace(-11, -4, 101), log=True),
    mixpods.plot_distributions(tree, "eps", bins=np.linspace(-11, -4, 101), log=True),
    mixpods.plot_distributions(
        tree, "ocean_vertical_heat_diffusivity", bins=np.linspace(-8, -1, 101), log=True
    ),
    # plot_distributions(tree, "Jq", bins=np.linspace(-1000, 0, 51), log=False),
    mixpods.plot_distributions(tree, "Rig_T", np.linspace(-0.5, 1.5, 61))
    * hv.VLine(0.25).opts(line_color="k"),
]
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:761: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:761: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
plot = (
    hv.Layout(handles)
    .opts("Overlay", frame_width=500)
    .cols(2)
    .opts(*mixpods.HV_TOOLS_OPTIONS)
)
select_legends(plot, figure_index=1, legend_position="right")
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:967: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(