This page was generated from docs/source/examples/example_2D_CROCO.ipynb.
2D Structure Function / CROCO Data
This example will guide you through each step necessary to compute scalar-based structure functions from a 2D CROCO data.
\(\textbf{General procedure}:\)
1 - Load a croco dataset, the dataset can be accessed https://doi.org/10.1175/JPO-D-18-0029.1
2 - Format dataset to match pyturbo_sf inputs
3 - Calculate second-order 2D longitudinal structure function as a function of separation distances
4 - Plot the 2D SF with errors and number of bootstraps
5 - Calculate Isotropic second-order longitudinal & Transverse structure function and their associated errors
6 - Plot Isotropic Longitudinal and transverse SFs with isotropy and homogeneity errors
7 - Calculate Isotropic advective-SF \(\delta u.\delta ADV\)
8 - Plot Isotropic advective-SF \(\delta u.\delta ADV\) with isotropy and homogeneity errors
8 - Calculate Isotropic advective vorticity-SF \(\delta \omega_{z}.\delta ADV_{\omega_{z}}\)
9 - Plot Isotropic advective-SF \(\delta \omega_{z}.\delta ADV_{\omega_{z}}\) with isotropy and homogeneity errors
10 - Calculate Isotropic scalar-scalar \(\delta u \delta T\)
11 - Plot Isotropic scalar-scalar \(\delta u \delta T\) with isotropy and homogeneity errors
\(\textbf{Note:}\) Pick A bootsize which divides the data size by power of 2
\(\textbf{Motivation:}\) This example provides an ideal comparison to the surface drifter structure function analysis presented in Pearson et al. (2019), who demonstrated significant biases between Lagrangian and Eulerian structure functions in the Gulf of Mexico. While their study relied on computationally expensive Lagrangian drifter deployments and the CROCO model data offers a cost-effective alternative for computing both traditional and advective structure functions. The advective structure functions implemented here follow the methodology that proved particularly valuable for diagnosing energy transfer in convergent regions where traditional structure functions become biased. This computational approach enables systematic exploration of submesoscale dynamics without the prohibitive costs of extensive field campaigns, while still capturing the essential physics of turbulent cascades in regions of strong surface convergence and divergence.
\(\textbf{Reference:}\)
Pearson, J., Fox-Kemper, B., Barkan, R., Choi, J., Bracco, A., & McWilliams, J. C. (2019). Impacts of convergence on structure functions from surface drifters in the Gulf of Mexico. Journal of Physical Oceanography, 49(3), 675-690. https://doi.org/10.1175/JPO-D-18-0029.1
[1]:
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm, SymLogNorm
import matplotlib.ticker as ticker
import pyturbo_sf as psf
import pyproj
import numpy as np
linewidth = 2
fontsize = 12
plt.rcParams['xtick.labelsize'] = fontsize
plt.rcParams['ytick.labelsize'] = fontsize
plt.rcParams['xtick.major.width'] = 2
plt.rcParams['xtick.minor.width'] = 2
plt.rcParams['ytick.major.width'] = 2
plt.rcParams['ytick.minor.width'] = 2
plt.rcParams['xtick.major.size'] = 10
plt.rcParams['xtick.minor.size'] = 5
plt.rcParams['ytick.major.size'] = 10
plt.rcParams['ytick.minor.size'] = 5
plt.rcParams['savefig.dpi'] = 150
plt.rc('font', family='serif')
import gc
Load the CROCO Dataset
[2]:
path_data = '/home/aayouche/Downloads/'
idt = 10
dsh = xr.open_dataset(path_data+'/'+'surface.nc')
dsg = xr.open_dataset(path_data+'/'+'GOM_150_grd.nc')
lon = dsg.lon_rho
lat = dsg.lat_rho
lon_center = np.nanmin(lon)
utm_zone = int(1 + (lon_center + 180.0) / 6.0)
# Create the projection
proj = pyproj.Proj(proj='utm', zone=utm_zone, ellps='WGS84', datum='WGS84')
x, y = proj(lon.values, lat.values)
X = 0.25*(x[1:,1:] + x[1:,:-1] + x[:-1,:-1] + x[:-1,1:])
Y = 0.25*(y[1:,1:] + y[1:,:-1] + y[:-1,:-1] + y[:-1,1:])
u = np.array(dsh.u.isel(s_rho=0,time=idt))
u = 0.5*(u[1:,:] + u[:-1,:])
v = np.array(dsh.v.isel(s_rho=0,time=idt))
v = 0.5*(v[:,1:] + v[:,:-1])
dudx = np.gradient(u,axis=1)/(2*150.)
dvdx = np.gradient(v,axis=1)/(2*150.)
dudy = np.gradient(u,axis=0)/(2*150.)
dvdy = np.gradient(v,axis=0)/(2*150.)
adv_u = u*dudx + v*dudy
adv_v = u*dvdx + v*dvdy
temp = np.array(dsh.temp.isel(s_rho=0,time=idt))
temp = 0.25*(temp[1:,1:] + temp[1:,:-1] + temp[:-1,:-1] + temp[:-1,1:])
f = 1.0e-4
vort = dvdx - dudy + f
dvortdx = np.gradient(vort,axis=1)/(2*150.)
dvortdy = np.gradient(vort,axis=0)/(2*150.)
adv_vort = u*dvortdx + v*dvortdy
gc.collect()
[2]:
70
[ ]:
Make A quick Plot
[3]:
fig,ax = plt.subplots(figsize=(6,6))
im = ax.pcolormesh(X/1.0e3,Y/1.0e3,0.5*(u**2 + v**2),cmap='jet',shading='gouraud')
pos1cb = ax.get_position()
cbar_ax = fig.add_axes([pos1cb.x0+pos1cb.width/4., pos1cb.y0+0.04+pos1cb.height , pos1cb.width/2. , 0.009])
cbar = fig.colorbar(im, cax=cbar_ax,orientation='horizontal',\
extend='both',ticks=np.arange(0,1,0.2))
cbar.set_label(r'$ {\rm KE \, [} {\rm m^{2}~} {\rm s^{-2}]}$',fontsize= fontsize + 2 ,labelpad=10)
cbar.ax.tick_params(labelsize=fontsize,direction='in',width=linewidth,size=4)
cbar.ax.xaxis.set_label_position('top')
cbar.ax.xaxis.set_ticks_position('both')
cbar.outline.set_linewidth( linewidth )
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth( linewidth )
ax.spines[axis].set_color('k')
ax.tick_params(direction='in',bottom='on',top='on',left='on',right='on')
ax.set_xlabel('X [km]')
ax.set_ylabel('Y [km]')
[3]:
Text(0, 0.5, 'Y [km]')
[4]:
gc.collect()
[4]:
3323
Format Dataset
[5]:
# Create the xarray dataset
ds = xr.Dataset({
'u': (['y', 'x'], u),
'v': (['y', 'x'], v),
'adv_u':(['y', 'x'], adv_u),
'adv_v':(['y', 'x'], adv_v),
'temperature':(['y', 'x'], temp),
'vort':(['y', 'x'], vort),
'adv_vort':(['y', 'x'], adv_vort),
}, coords={
'x': (['y', 'x'], X),
'y': (['y', 'x'], Y)
})
[6]:
ds
[6]:
<xarray.Dataset> Size: 69MB
Dimensions: (y: 1537, x: 1025)
Coordinates:
y (y, x) float64 13MB 3.039e+06 3.039e+06 ... 3.267e+06 3.267e+06
x (y, x) float64 13MB 3.542e+05 3.544e+05 ... 5.086e+05 5.088e+05
Data variables:
u (y, x) float32 6MB -0.3907 -0.4025 -0.4012 ... 0.01855 0.01787
v (y, x) float32 6MB 0.08256 0.09225 0.1009 ... 0.0189 0.01809
adv_u (y, x) float32 6MB 1.877e-05 1.042e-05 ... 2.56e-07 9.305e-08
adv_v (y, x) float32 6MB -1.566e-05 -1.693e-05 ... -3.382e-07
temperature (y, x) float32 6MB 25.06 25.06 25.06 ... 20.97 20.97 20.97
vort (y, x) float32 6MB 9.174e-05 9.375e-05 ... 7.981e-05 8.993e-05
adv_vort (y, x) float32 6MB -6.696e-09 -4.602e-09 ... 7.135e-10Keywords for 2D and Isotropic SF
[7]:
# 2D Binning
sf_kwargs_2d = dict(
ds=ds,
variables_names=["u", "v"],
order=2,
bins={
'x': np.logspace(np.log10(1.0e3), np.log10(150.0e3), 12),
'y': np.logspace(np.log10(1.0e3), np.log10(150.0e3), 12)
},
bootsize={'y': 96, 'x': 64},
initial_nbootstrap=60,
max_nbootstrap=420,
step_nbootstrap=180,
n_jobs=-1,
backend='loky',
seed=42
)
# Isotropic
sf_kwargs_is = dict(
ds=ds,
bins={
'r': np.logspace(np.log10(1.0e3), np.log10(150.0e3), 21)
},
bootsize={'y': 12, 'x': 8},
initial_nbootstrap=100000,
max_nbootstrap=100000,
step_nbootstrap=0,
confidence_interval=0.95,
n_jobs=-1,
backend='loky',
seed=42
)
Calculate 2D Longitudinal SF, \(\epsilon=0.05\)
[8]:
%%time
sf_result = psf.bin_sf_2d(**sf_kwargs_2d, fun = 'longitudinal',convergence_eps=0.05)
Dimensions ('y', 'x') are already in the expected order
Using bootsize: {'y': 96, 'x': 64}
Bootstrappable dimensions: ['y', 'x']
Two bootstrappable dimensions. Available spacings: [1, 2, 4, 8, 16]
============================================================
STARTING BIN_SF WITH FUNCTION TYPE: longitudinal
Variables: ['u', 'v'], Order: 2
Using seed 42 for reproducible bootstrap sampling
============================================================
INITIAL BOOTSTRAP PHASE
Processing spacing 1 with 12 bootstraps
Processing spacing 2 with 12 bootstraps
Processing spacing 4 with 12 bootstraps
Processing spacing 8 with 12 bootstraps
Processing spacing 16 with 12 bootstraps
CALCULATING BIN DENSITIES
Total points collected: 580700160
Bins with points: 121/121
Marked 118 bins as converged (converged_eps)
STARTING ADAPTIVE CONVERGENCE LOOP
Iteration 1 - 3 unconverged bins
Grouped unconverged bins into 1 groups
Processing 3 bins with step size 180 in density quartile 0
Bin (8,1) CONVERGED with std 0.007322
Bin (9,1) CONVERGED with std 0.005053
Bin (10,1) CONVERGED with std 0.008340
All bins have converged or reached max bootstraps!
FINAL CONVERGENCE STATISTICS:
Total bins with data (>10 points): 121
Converged bins: 121
Unconverged bins: 0
Bins at max bootstraps: 0
Creating output dataset...
2D SF COMPLETED SUCCESSFULLY!
============================================================
CPU times: user 2.52 s, sys: 448 ms, total: 2.97 s
Wall time: 19.2 s
[9]:
gc.collect()
[9]:
37
Plot 2D Longitudial SF (structure function, standard error, number of bootstrap samples to reach convergence)
[10]:
# Create figure with 3 subplots
fig, axes = plt.subplots(1, 3, figsize=(16, 5.5), constrained_layout=True)
# Common settings for all plots
plot_titles = ['Structure Function', 'Standard Error', 'Bootstraps Required']
cmaps = ['viridis', 'viridis', 'plasma']
labels = [r'$\delta u_{LL}$ [$m^{2}~s^{-2}$]', r'$\delta u_{LL}$ [$m^{2}~s^{-2}$]', 'Number of Bootstraps']
# Data to plot
data_to_plot = [sf_result.sf, sf_result.std_error, sf_result.nbootstraps]
for i, (ax, data, title, cmap, label) in enumerate(zip(axes, data_to_plot, plot_titles, cmaps, labels)):
# Create the plot with specific shading for the first plot
if i == 0 or i == 1:
im = ax.pcolormesh(sf_result.x, sf_result.y, data, cmap=cmap, shading='gouraud')
else:
im = ax.pcolormesh(sf_result.x, sf_result.y, data, cmap=cmap)
# Add colorbar with proper formatting
cbar = fig.colorbar(im, ax=ax)
cbar.set_label(label, fontsize=11, fontweight='bold')
# Set axis labels and title
ax.set_xlabel(r'$r_x$ [m]', fontsize=11, fontweight='bold')
if i == 0:
ax.set_ylabel(r'$r_y$ [m]', fontsize=11, fontweight='bold')
# Set log scales
ax.set_xscale('log')
ax.set_yscale('log')
# Add title
ax.set_title(title, fontsize=12, fontweight='bold', pad=10)
# Customize spines
for spine in ax.spines.values():
spine.set_linewidth(1.5)
# Customize tick parameters
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
# Add scientific notation for axes
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
ax.yaxis.set_major_formatter(ticker.LogFormatterSciNotation())
# Add main title for the entire figure
fig.suptitle('CROCO Structure Function Analysis', fontsize=14, fontweight='bold', y=1.08)
# Add annotation explaining the analysis
fig.text(0.5, -0.08,
'Structure function analysis reveals spatial variability of the longitudinal SF across scales.\n'
'Threshold value of 0.05 $m^{2}~s^{-2}$ highlights Gulf Stream frontal features.',
ha='center', fontsize=10, fontstyle='italic')
plt.show()
Physical Interpretation: 2D Structure Functions
The 2D structure function maps reveal the anisotropic structure of ocean turbulence:
Structure Function panel: Shows how velocity correlations vary with both separation distance and direction. Anisotropy appears as elongation along preferred axes (e.g., along-front vs cross-front).
Standard Error panel: Quantifies uncertainty in each bin. Larger errors at the domain edges reflect fewer sample pairs.
Bootstrap Iterations panel: Shows where more computational effort was needed to achieve convergence, typically at scales with high intermittency or sparse sampling.
[11]:
gc.collect()
[11]:
2561
Calculate ISOTROPIC Longitudinal and Transverse SF, \(\epsilon=0.05\)
[12]:
%%time
iso_sf_lon = psf.get_isotropic_sf_2d(**sf_kwargs_is, variables_names=["u","v"], order=2, fun='longitudinal', convergence_eps=0.05)
iso_sf_trans = psf.get_isotropic_sf_2d(**sf_kwargs_is, variables_names=["u","v"], order=2, fun='transverse', convergence_eps=0.05)
Dimensions ('y', 'x') are already in the expected order
Using bootsize: {'y': 12, 'x': 8}
Bootstrappable dimensions: ['y', 'x']
Two bootstrappable dimensions. Available spacings: [1, 2, 4, 8, 16, 32, 64, 128]
============================================================
STARTING ISOTROPIC_SF WITH FUNCTION TYPE: longitudinal
Variables: ['u', 'v'], Order: 2
Confidence level: 0.95
============================================================
INITIAL BOOTSTRAP PHASE
Processing spacing 1 with 12500 bootstraps
Processing spacing 2 with 12500 bootstraps
Processing spacing 4 with 12500 bootstraps
Processing spacing 8 with 12500 bootstraps
Processing spacing 16 with 12500 bootstraps
Processing spacing 32 with 12500 bootstraps
Processing spacing 64 with 12500 bootstraps
Processing spacing 128 with 12500 bootstraps
CALCULATING BIN DENSITIES
Total points collected: 271200000
Bins with points: 20/20
Marked 20 bins as converged (converged_eps)
STARTING ADAPTIVE CONVERGENCE LOOP
All bins have converged or reached max bootstraps!
FINAL CONVERGENCE STATISTICS:
Total bins with data (>10 points): 20
Converged bins: 20
Unconverged bins: 0
Bins at max bootstraps: 20
Creating output dataset...
ISOTROPIC SF COMPLETED SUCCESSFULLY!
============================================================
Dimensions ('y', 'x') are already in the expected order
Using bootsize: {'y': 12, 'x': 8}
Bootstrappable dimensions: ['y', 'x']
Two bootstrappable dimensions. Available spacings: [1, 2, 4, 8, 16, 32, 64, 128]
============================================================
STARTING ISOTROPIC_SF WITH FUNCTION TYPE: transverse
Variables: ['u', 'v'], Order: 2
Confidence level: 0.95
============================================================
INITIAL BOOTSTRAP PHASE
Processing spacing 1 with 12500 bootstraps
Processing spacing 2 with 12500 bootstraps
Processing spacing 4 with 12500 bootstraps
Processing spacing 8 with 12500 bootstraps
Processing spacing 16 with 12500 bootstraps
Processing spacing 32 with 12500 bootstraps
Processing spacing 64 with 12500 bootstraps
Processing spacing 128 with 12500 bootstraps
CALCULATING BIN DENSITIES
Total points collected: 271200000
Bins with points: 20/20
Marked 20 bins as converged (converged_eps)
STARTING ADAPTIVE CONVERGENCE LOOP
All bins have converged or reached max bootstraps!
FINAL CONVERGENCE STATISTICS:
Total bins with data (>10 points): 20
Converged bins: 20
Unconverged bins: 0
Bins at max bootstraps: 20
Creating output dataset...
ISOTROPIC SF COMPLETED SUCCESSFULLY!
============================================================
CPU times: user 1min 2s, sys: 3.06 s, total: 1min 6s
Wall time: 1min 34s
Plot Isotropic Longitudinal and Transverse SF
[13]:
# Create separate figures for better spacing and clarity
plt.figure(figsize=(10, 6))
# Structure Function with Confidence Interval
# Add longitudinal structure function with error bars
plt.loglog(iso_sf_lon.r, iso_sf_lon.sf, 'b-', linewidth=2.5, label='Longitudinal Structure Function'+' '+'$\delta u_{LL}$')
plt.fill_between(iso_sf_lon.r,iso_sf_lon.ci_lower,iso_sf_lon.ci_upper,color='b',alpha=0.3)
# Add transverse structure function with error bars
plt.loglog(iso_sf_trans.r, iso_sf_trans.sf, 'r-', linewidth=2.5, label='Transverse Structure Function'+' '+'$\delta u_{TT}$')
plt.fill_between(iso_sf_trans.r,iso_sf_trans.ci_lower,iso_sf_trans.ci_upper,color='r',alpha=0.3)
# Formatting
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel(r'$[m^{2}~s^{-2}]$', fontsize=12, fontweight='bold')
plt.title('Isotropic Longitudinal & Transverse Structure Functions with Confidence Intervals',
fontsize=14, fontweight='bold', pad=10)
plt.legend(fontsize=11, frameon=True, facecolor='white', framealpha=0.9)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.yscale('log')
plt.xscale('log')
# Customize spines
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
# Customize tick parameters
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
# --------------------------------------------------
# Errors of Isotropy and Homogeneity
# --------------------------------------------------
plt.figure(figsize=(12, 5))
# Error of isotropy
plt.subplot(121)
plt.plot(iso_sf_lon.r, iso_sf_lon.error_isotropy, 'b-', linewidth=2.5,label='Longitudinal Structure Function'+' '+'$\delta u_{LL}$')
plt.plot(iso_sf_trans.r, iso_sf_trans.error_isotropy, 'r-', linewidth=2.5,label='Transverse Structure Function'+' '+'$\delta u_{TT}$')
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Isotropy'+' '+r'$[m^{2}~s^{-2}]$', fontsize=12, fontweight='bold')
plt.title('Error of Isotropy (Variation with Angle)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
plt.legend(fontsize=11, frameon=True, facecolor='white', framealpha=0.9)
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
# Error of homogeneity
plt.subplot(122)
plt.plot(iso_sf_lon.r_subset, iso_sf_lon.error_homogeneity, 'b-', linewidth=2.5,label='Longitudinal Structure Function'+' '+'$\delta u_{LL}$')
plt.plot(iso_sf_trans.r_subset, iso_sf_trans.error_homogeneity, 'r-', linewidth=2.5,label='Transverse Structure Function'+' '+'$\delta u_{TT}$')
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Homogeneity'+' '+r'$[m^{2}~s^{-2}]$', fontsize=12, fontweight='bold')
plt.title('Error of Homogeneity (Variation with Radius)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
plt.legend(fontsize=11, frameon=True, facecolor='white', framealpha=0.9)
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
plt.show()
Physical Interpretation: Isotropic Structure Functions
The isotropic longitudinal and transverse structure functions characterize energy distribution across scales:
Longitudinal SF (\(S_{LL}\)): Measures velocity variance in the direction of separation. In 2D turbulence, the ratio of the slope to theoretical predictions (-2/3 for enstrophy cascade, -2 for energy cascade) identifies the dominant cascade regime.
Transverse SF (\(S_{TT}\)): Measures velocity variance perpendicular to separation. For incompressible 2D flow, \(S_{TT}/S_{LL} \approx 2\) indicates isotropy.
Deviations from this ratio reveal anisotropic forcing or boundary effects.
Calculate Isotropic advective-SF \(\delta u.\delta ADV\)
[14]:
%%time
iso_sf_adv = psf.get_isotropic_sf_2d(**sf_kwargs_is,
order=1,
variables_names=["u","v","adv_u","adv_v"],
fun='advective',
convergence_eps=1.0e-8,
)
Dimensions ('y', 'x') are already in the expected order
Using bootsize: {'y': 12, 'x': 8}
Bootstrappable dimensions: ['y', 'x']
Two bootstrappable dimensions. Available spacings: [1, 2, 4, 8, 16, 32, 64, 128]
============================================================
STARTING ISOTROPIC_SF WITH FUNCTION TYPE: advective
Variables: ['u', 'v', 'adv_u', 'adv_v'], Order: 1
Confidence level: 0.95
============================================================
INITIAL BOOTSTRAP PHASE
Processing spacing 1 with 12500 bootstraps
Processing spacing 2 with 12500 bootstraps
Processing spacing 4 with 12500 bootstraps
Processing spacing 8 with 12500 bootstraps
Processing spacing 16 with 12500 bootstraps
Processing spacing 32 with 12500 bootstraps
Processing spacing 64 with 12500 bootstraps
Processing spacing 128 with 12500 bootstraps
CALCULATING BIN DENSITIES
Total points collected: 271200000
Bins with points: 20/20
Marked 20 bins as converged (converged_eps)
STARTING ADAPTIVE CONVERGENCE LOOP
All bins have converged or reached max bootstraps!
FINAL CONVERGENCE STATISTICS:
Total bins with data (>10 points): 20
Converged bins: 20
Unconverged bins: 0
Bins at max bootstraps: 20
Creating output dataset...
ISOTROPIC SF COMPLETED SUCCESSFULLY!
============================================================
CPU times: user 31.5 s, sys: 1.42 s, total: 32.9 s
Wall time: 49.4 s
Plot Isotropic advective-SF \(\delta u.\delta ADV\)
[15]:
%%time
# Create separate figures for better spacing and clarity
plt.figure(figsize=(10, 6))
# Structure Function with Confidence Interval
plt.plot(iso_sf_adv.r, iso_sf_adv.sf, 'b-', linewidth=2.5)
plt.fill_between(iso_sf_adv.r,
iso_sf_adv.ci_lower,
iso_sf_adv.ci_upper,
alpha=0.3, color='blue')
# Formatting
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel(r'$\delta u.\delta ADV~[m^{2}~s^{-3}]$', fontsize=12, fontweight='bold')
plt.title('Isotropic Advective Structure Function with Confidence Intervals',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
# Customize tick parameters
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
# --------------------------------------------------
# Errors of Isotropy and Homogeneity
# --------------------------------------------------
plt.figure(figsize=(12, 5))
# Error of isotropy
plt.subplot(121)
plt.plot(iso_sf_adv.r, iso_sf_adv.error_isotropy, 'b-', linewidth=2.5)
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Isotropy'+' '+r'$[m^{2}~s^{-3}]$', fontsize=12, fontweight='bold')
plt.title('Error of Isotropy (Variation with Angle)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
# Error of homogeneity
plt.subplot(122)
plt.plot(iso_sf_adv.r_subset, iso_sf_adv.error_homogeneity, 'b-', linewidth=2.5)
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Homogeneity'+' '+r'$[m^{2}~s^{-3}]$', fontsize=12, fontweight='bold')
plt.title('Error of Homogeneity (Variation with Radius)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
plt.show()
CPU times: user 378 ms, sys: 8.9 ms, total: 387 ms
Wall time: 384 ms
Physical Interpretation: Velocity Advective Structure Function
The advective SF \(\langle \delta \mathbf{u} \cdot \delta \mathbf{A}_u \rangle\) directly measures kinetic energy flux:
Negative values indicate inverse cascade (energy moving to larger scales), characteristic of quasi-geostrophic turbulence.
Scale-independent plateau in the inertial range confirms constant energy flux.
This diagnostic is less sensitive to anisotropy than \(S_{LLL}\), making it more robust for real ocean data.
Calculate Isotropic \(\delta \omega_{z} * \delta ADV_{\omega_{z}}\)
[16]:
%%time
iso_sf_adv = psf.get_isotropic_sf_2d(**sf_kwargs_is,
variables_names=["vort","adv_vort"],
order=(1,1),
fun='scalar_scalar',
convergence_eps=1.0e-14
)
Dimensions ('y', 'x') are already in the expected order
Using bootsize: {'y': 12, 'x': 8}
Bootstrappable dimensions: ['y', 'x']
Two bootstrappable dimensions. Available spacings: [1, 2, 4, 8, 16, 32, 64, 128]
============================================================
STARTING ISOTROPIC_SF WITH FUNCTION TYPE: scalar_scalar
Variables: ['vort', 'adv_vort'], Order: (1, 1)
Confidence level: 0.95
============================================================
INITIAL BOOTSTRAP PHASE
Processing spacing 1 with 12500 bootstraps
Processing spacing 2 with 12500 bootstraps
Processing spacing 4 with 12500 bootstraps
Processing spacing 8 with 12500 bootstraps
Processing spacing 16 with 12500 bootstraps
Processing spacing 32 with 12500 bootstraps
Processing spacing 64 with 12500 bootstraps
Processing spacing 128 with 12500 bootstraps
CALCULATING BIN DENSITIES
Total points collected: 271200000
Bins with points: 20/20
Marked 20 bins as converged (converged_eps)
STARTING ADAPTIVE CONVERGENCE LOOP
All bins have converged or reached max bootstraps!
FINAL CONVERGENCE STATISTICS:
Total bins with data (>10 points): 20
Converged bins: 20
Unconverged bins: 0
Bins at max bootstraps: 20
Creating output dataset...
ISOTROPIC SF COMPLETED SUCCESSFULLY!
============================================================
CPU times: user 31.6 s, sys: 1.42 s, total: 33.1 s
Wall time: 43.8 s
Plot Isotropic \(\delta \omega_{z} * \delta ADV_{\omega_{z}}\)
[17]:
# Create separate figures for better spacing and clarity
plt.figure(figsize=(10, 6))
# Structure Function with Confidence Interval
plt.plot(iso_sf_adv.r, iso_sf_adv.sf, 'b-', linewidth=2.5)
plt.fill_between(iso_sf_adv.r,
iso_sf_adv.ci_lower,
iso_sf_adv.ci_upper,
alpha=0.3, color='blue')
# Formatting
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel(r'$\delta \omega_{z}\delta ADV_{\omega}~[s^{-3}]$', fontsize=12, fontweight='bold')
plt.title('Isotropic Advective-Vorticity Structure Function with Confidence Intervals',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
# Customize tick parameters
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
# --------------------------------------------------
# Errors of Isotropy and Homogeneity
# --------------------------------------------------
plt.figure(figsize=(12, 5))
# Error of isotropy
plt.subplot(121)
plt.plot(iso_sf_adv.r, iso_sf_adv.error_isotropy, 'b-', linewidth=2.5)
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Isotropy'+' '+r'$[s^{-3}]$', fontsize=12, fontweight='bold')
plt.title('Error of Isotropy (Variation with Angle)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
# Error of homogeneity
plt.subplot(122)
plt.plot(iso_sf_adv.r_subset, iso_sf_adv.error_homogeneity, 'b-', linewidth=2.5)
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Homogeneity'+' '+r'$[s^{-3}]$', fontsize=12, fontweight='bold')
plt.title('Error of Homogeneity (Variation with Radius)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
plt.show()
Physical Interpretation: Vorticity Advective Structure Function
The vorticity-based advective SF \(\langle \delta \omega \cdot \delta A_\omega \rangle\) measures enstrophy (squared vorticity) flux:
In 2D turbulence, enstrophy cascades forward to smaller scales while energy cascades inversely to larger scales.
Positive values here indicate forward enstrophy cascade, complementing the inverse energy cascade seen in velocity SFs.
The enstrophy flux rate \(\eta\) relates to this SF, providing the other half of the dual cascade picture.
Calculate Isotropic \(\delta u \delta T\)
[18]:
%%time
iso_sf_scalar = psf.get_isotropic_sf_2d(**sf_kwargs_is,
variables_names=["u","temperature"],
order=(1,1),
fun='scalar_scalar',
convergence_eps=0.01,
)
Dimensions ('y', 'x') are already in the expected order
Using bootsize: {'y': 12, 'x': 8}
Bootstrappable dimensions: ['y', 'x']
Two bootstrappable dimensions. Available spacings: [1, 2, 4, 8, 16, 32, 64, 128]
============================================================
STARTING ISOTROPIC_SF WITH FUNCTION TYPE: scalar_scalar
Variables: ['u', 'temperature'], Order: (1, 1)
Confidence level: 0.95
============================================================
INITIAL BOOTSTRAP PHASE
Processing spacing 1 with 12500 bootstraps
Processing spacing 2 with 12500 bootstraps
Processing spacing 4 with 12500 bootstraps
Processing spacing 8 with 12500 bootstraps
Processing spacing 16 with 12500 bootstraps
Processing spacing 32 with 12500 bootstraps
Processing spacing 64 with 12500 bootstraps
Processing spacing 128 with 12500 bootstraps
CALCULATING BIN DENSITIES
Total points collected: 271200000
Bins with points: 20/20
Marked 20 bins as converged (converged_eps)
STARTING ADAPTIVE CONVERGENCE LOOP
All bins have converged or reached max bootstraps!
FINAL CONVERGENCE STATISTICS:
Total bins with data (>10 points): 20
Converged bins: 20
Unconverged bins: 0
Bins at max bootstraps: 20
Creating output dataset...
ISOTROPIC SF COMPLETED SUCCESSFULLY!
============================================================
CPU times: user 31.5 s, sys: 1.49 s, total: 33 s
Wall time: 43.7 s
Plot Isotropic \(\delta u \delta T\)
[19]:
# Create separate figures for better spacing and clarity
plt.figure(figsize=(10, 6))
# Structure Function with Confidence Interval
plt.plot(iso_sf_scalar.r, iso_sf_scalar.sf, 'b-', linewidth=2.5)
plt.fill_between(iso_sf_scalar.r,
iso_sf_scalar.ci_lower,
iso_sf_scalar.ci_upper,
alpha=0.3, color='blue')
# Formatting
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel(r'$\delta u \delta T~[m~s^{-1}~°C]$', fontsize=12, fontweight='bold')
plt.title('Isotropic Velocity-Scalar Structure Function with Confidence Intervals',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
# Customize tick parameters
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
# --------------------------------------------------
# Errors of Isotropy and Homogeneity
# --------------------------------------------------
plt.figure(figsize=(12, 5))
# Error of isotropy
plt.subplot(121)
plt.plot(iso_sf_scalar.r, iso_sf_scalar.error_isotropy, 'b-', linewidth=2.5)
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Isotropy'+' '+r'$[m~s^{-1}~°C]$', fontsize=12, fontweight='bold')
plt.title('Error of Isotropy (Variation with Angle)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
# Error of homogeneity
plt.subplot(122)
plt.plot(iso_sf_scalar.r_subset, iso_sf_scalar.error_homogeneity, 'b-', linewidth=2.5)
plt.xlabel(r'$r$ [m]', fontsize=12, fontweight='bold')
plt.ylabel('Error of Homogeneity'+' '+r'$[m~s^{-1}~°C]$', fontsize=12, fontweight='bold')
plt.title('Error of Homogeneity (Variation with Radius)',
fontsize=14, fontweight='bold', pad=10)
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.xscale('log')
# Customize spines and ticks
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(1.5)
ax.tick_params(which='both', direction='in', width=1.2, length=5,
labelsize=10, top=True, right=True)
ax.xaxis.set_major_formatter(ticker.LogFormatterSciNotation())
plt.tight_layout()
plt.show()
Physical Interpretation: Velocity-Temperature Cross Structure Function
The cross-SF \(\langle \delta u \cdot \delta T \rangle\) quantifies the correlation between velocity and temperature fluctuations:
Non-zero values indicate turbulent heat transport — velocity fluctuations systematically correlated with temperature anomalies.
The sign indicates the direction of heat flux: positive for downgradient transport (warm water moving with the flow).
Scale dependence reveals which eddy sizes contribute most to cross-front heat exchange.
This is particularly relevant for understanding mesoscale eddy heat transport in the ocean.
[ ]:
[ ]: