Contents

Heat Transport

We will assume the heat budget in the model takes the form

(32)\[\begin{align} \frac{D \Theta}{D t} = - \nabla \cdot (\mathbf{F}^\Theta_{diff}) - \frac{1}{\rho_0 c_p }\frac{\partial Q}{\partial z} \end{align}\]

The term on the LHS represents the advective flux of conservative temperature. The first term on the RHS is the convergence of diffusive fluxes of heat. In the model, this represents turbulent diffusion due to mesoscale eddies and small-scale mixing. The final term is the convergence of heat fluxes from the atmosphere.

We will further assume a rigid lid. Note that the ECCO model actually uses a sophisticated nonlinear free surface form of the conservation equations, so these approximations aren’t quite right. But our goal here is to focus on the big picture, and these details would be a distraction.

import xgcm
import xarray as xr
from matplotlib import pyplot as plt
import numpy as np
import hvplot.xarray

from ecco_utils import load_ecco_v4

plt.rcParams['figure.figsize'] = (12, 6)

ds = load_ecco_v4()
ds

<xarray.Dataset>
Dimensions:    (face: 13, lon_c: 360, lon_g: 360, lat_c: 240, lat_g: 240, Z: 50, Zl: 50, Zp1: 51, Zu: 50, time: 288, time_snp: 287)
Coordinates: (12/42)
  * face       (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
    i          (lon_c) int64 0 1 2 3 4 5 6 7 ... 352 353 354 355 356 357 358 359
    i_g        (lon_g) int64 0 1 2 3 4 5 6 7 ... 352 353 354 355 356 357 358 359
    j          (lat_c) int64 30 31 32 33 34 35 36 ... 264 265 266 267 268 269
    j_g        (lat_g) int64 30 31 32 33 34 35 36 ... 264 265 266 267 268 269
    k          (Z) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    dyG        (lat_c, lon_g) float32 dask.array<chunksize=(60, 90), meta=np.ndarray>
    dxG        (lat_g, lon_c) float32 dask.array<chunksize=(60, 90), meta=np.ndarray>
  * lon_c      (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lon_g      (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c      (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lat_g      (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables: (12/35)
    ADVr_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    ADVr_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrE_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrE_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrI_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrI_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    ...         ...
    UVELMASS   (time, Z, lat_c, lon_g) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    VVELMASS   (time, Z, lat_g, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    UVELSTAR   (time, Z, lat_c, lon_g) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    VVELSTAR   (time, Z, lat_g, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    oceTAUX    (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    oceTAUY    (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>

Vertically Integrated Heat Budget

Integrating the budget vertically and assuming \(\eta = 0\), we obtain

(33)\[\begin{align} \frac{\partial}{\partial t} \int_{-H}^0 \Theta dz = - \frac{\partial}{\partial x} \int_{-H}^0 u \Theta dz - \frac{\partial}{\partial y} \int_{-H}^0 v \Theta dz - \frac{\partial}{\partial x} \int_{-H}^0 F^\Theta_{diff,x} dz - \frac{\partial}{\partial y} \int_{-H}^0 F^\Theta_{diff,y} dz + \frac{Q}{\rho_0 c_p} \end{align}\]

We can write this in a more compact form as

(34)\[\begin{align} \frac{\partial}{\partial t} \int_{-H}^0 \Theta dz = - \nabla_H \cdot (\mathbf{F}^\Theta_{adv,H} + \mathbf{F}^\Theta_{diff,H}) + \frac{Q}{\rho_0 c_p} \end{align}\]

Where \(\mathbf{F}^\Theta_{adv,H}\) and \(\mathbf{F}^\Theta_{diff,H}\) are the vertically integrated advective and diffusive (2D) heat flux vectors. The operator \(\nabla_H\) represents the divergence in two dimensions (x, y).

from dask_gateway import Gateway
gway = Gateway()
cluster = gway.new_cluster()
cluster

/srv/conda/envs/notebook/lib/python3.8/site-packages/dask_gateway/client.py:21: FutureWarning: format_bytes is deprecated and will be removed in a future release. Please use dask.utils.format_bytes instead.
  from distributed.utils import LoopRunner, format_bytes

grid = xgcm.Grid(ds, periodic=['X'])

fluxes = xr.Dataset({
    'F_adv_x': ds.ADVx_TH.sum('Z'),
    'F_adv_y': ds.ADVy_TH.sum('Z'),
    'F_diff_x': ds.DFxE_TH.sum('Z'),
    'F_diff_y': ds.DFyE_TH.sum('Z'),    
}).reset_coords(drop=True)
fluxes

<xarray.Dataset>
Dimensions:   (lon_g: 360, lat_c: 240, time: 288, lon_c: 360, lat_g: 240)
Coordinates:
  * time      (time) datetime64[ns] 1992-01-15 1992-02-13 ... 2015-12-14
  * lon_g     (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c     (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lon_c     (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lat_g     (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables:
    F_adv_x   (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    F_adv_y   (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    F_diff_x  (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    F_diff_y  (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
try:
    fluxes = xr.open_dataset('tmp_data/ECCO_lateral_heat_fluxes.nc')
except FileNotFoundError:
    import dask
    cluster.scale(20)
    with cluster.get_client():
        for vname in fluxes.data_vars:
            fluxes[vname].load()
    cluster.scale(0)
    fluxes.to_netcdf('tmp_data/ECCO_lateral_heat_fluxes.nc')
fluxes

<xarray.Dataset>
Dimensions:   (time: 288, lon_g: 360, lat_c: 240, lon_c: 360, lat_g: 240)
Coordinates:
  * time      (time) datetime64[ns] 1992-01-15 1992-02-13 ... 2015-12-14
  * lon_g     (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c     (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lon_c     (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lat_g     (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables:
    F_adv_x   (time, lat_c, lon_g) float32 ...
    F_adv_y   (time, lat_g, lon_c) float32 ...
    F_diff_x  (time, lat_c, lon_g) float32 ...
    F_diff_y  (time, lat_g, lon_c) float32 ...
over_rA = (ds.rA.compute()**-1).fillna(0.)

advective_flux_convergence = -(
    grid.diff(fluxes.F_adv_x, 'X', boundary='extend') +
    grid.diff(fluxes.F_adv_y, 'Y', boundary='extend')
) * over_rA

diffusive_flux_convergence = -(
    grid.diff(fluxes.F_diff_x, 'X', boundary='extend') +
    grid.diff(fluxes.F_diff_y, 'Y', boundary='extend')
) * over_rA

rho0 = 1029.
cp = 3850.

land_mask = (ds.hFacC[0] > 0).compute()
surface_flux = ds.TFLUX.compute().where(land_mask) / (rho0 * cp)

clim = (-1e-4, 1e-4)
big_plot = (
 advective_flux_convergence.hvplot.quadmesh('lon_c', 'lat_c', clim=clim, cmap='RdBu_r', label='Advective_Flux_Convergence') +
 diffusive_flux_convergence.hvplot.quadmesh('lon_c', 'lat_c', clim=clim, cmap='RdBu_r', label='Meridional_Flux_Convergence') +
 surface_flux.hvplot.quadmesh('lon_c', 'lat_c',clim=clim, cmap='RdBu_r', label='Surface_Flux') 
).cols(1)
big_plot