Global Heat Budget Closure

Contributors: Jan-Erik Tesdal, Ryan Abernathey and Ian Fenty

A major part of this tutorial is based on “A Note on Practical Evaluation of Budgets in ECCO Version 4 Release 3” by Christopher G. Piecuch (https://ecco.jpl.nasa.gov/drive/files/Version4/Release3/doc/v4r3_budgets_howto.pdf). Calculation steps and Python code presented here are converted from the MATLAB code presented in the above reference.

Objectives

Evaluating and closing the heat budget over the global ocean.

Introduction

The ocean heat content (OHC) variability is described here with potential temperature (\theta) which is given by the ECCOv4 diagnostic output THETA. The budget equation describing the change in \theta is evaluated in general as

\frac{\partial \theta}{\partial t} = -\nabla \cdot (\theta \mathbf{u})-\nabla\cdot\mathbf{F}_\textrm{diff}^{\theta}+{F}_\textrm{forc}^{\theta}

The heat budget includes the change in temperature over time (\frac{\partial \theta}{\partial t}), the convergence of heat advection (-\nabla \cdot (\theta \mathbf{u})) and heat diffusion (-\nabla\cdot\mathbf{F}_\textrm{diff}), plus downward heat flux from the atmosphere ({F}_\textrm{forc}). Note that in our definition {F}_\textrm{forc} contains both latent and sensible air-sea heat fluxes, longwave and shortwave radiation, as well as geothermal heat flux.

In the special case of ECCOv4, the heat budget is formulated as

\underbrace{\frac{\partial(s^*\theta)}{\partial t}}_{G^{\theta}_\textrm{total}} = \underbrace{-\nabla_{z^{*}} \cdot(s^*\theta\,\mathbf{v}_{res}) - \frac{\partial(\theta\,w_{res})}{\partial z^{*}}}_{G^{\theta}_\textrm{advection}}\underbrace{- s^* ({\nabla\cdot\mathbf{F}_\textrm{diff}^{\theta}})}_{G^{\theta}_\textrm{diffusion}} + \underbrace{s^* {F}_\textrm{forc}^{\theta}}_{G^{\theta}_\textrm{forcing}}

where z^{*} = \frac{z - \eta}{H + \eta}H and \nabla_{z^{*}}/\frac{\partial}{\partial z^{*}} are horizontal/vertical divergences in the z^* frame. Also note that the advection is now separated into horizontal (\mathbf{v}_{res}) and vertical (w_{res}) components, and there is a scaling factor (s^* = 1+ \frac{\eta}{H}) applied to the horizontal advection as well as the diffusion term (G^{\theta}_\textrm{diffusion}) and forcing term (G^{\theta}_\textrm{forcing}). s^* is a function of \eta which is the displacement of the ocean surface from its resting position of z=0 (i.e., sea height anomaly). H is the ocean depth. s^{*} comes from the coordinate transformation from z to z^* (Campin and Adcroft, 2004; Campin et al., 2004). See ECCOv4 Global Volume Budget Closure for a more detailed explanation of the z^* coordinate system.

Note that the velocity terms in the ECCOv4 heat budget equation (\mathbf{v}_{res} and w_{res}) are described as the “residual mean” velocities, which contain both the resolved (Eulerian) flow field, as well as the “GM bolus” velocity (i.e., parameterizing unresolved eddy effects):

(u_{res},v_{res},w_{res})= (u,v,w)+ (u_b,v_b,w_b)

Here (u_b,v_b,w_b) is the bolus velocity parameter, taking into account the correlation between velocity and thickness (also known as the eddy induced transportor the eddy advection term).

Evaluating the heat budget

We will evalute each term in the above heat budget

G^{\theta}_\textrm{total} = G^{\theta}_\textrm{advection} + G^{\theta}_\textrm{diffusion} + G^{\theta}_\textrm{forcing}

The total tendency of \theta (G^{\theta}_\textrm{total}) is the sum of the \theta tendencies from advective heat convergence (G^{\theta}_\textrm{advection}), diffusive heat convergence (G^{\theta}_\textrm{diffusion}) and total forcing (G^{\theta}_\textrm{forcing}).

We present calculation sequentially for each term starting with G^{\theta}_\textrm{total} which will be derived by differencing instantaneous monthly snapshots of \theta. The terms on the right hand side of the heat budget are derived from monthly-averaged fields.

Datasets to download

If you don’t have any of the following datasets already, you will need to download them to complete the tutorial. Aside from the grid geometry file (which has no time dimension), you will need 2 monthly datasets for the full time span of ECCOv4r4 output (1992 through 2017). The ShortNames of the datasets are:

  • ECCO_L4_GEOMETRY_LLC0090GRID_V4R4

  • ECCO_L4_OCEAN_3D_TEMPERATURE_FLUX_LLC0090GRID_MONTHLY_V4R4 (1993-2016)

  • ECCO_L4_HEAT_FLUX_LLC0090GRID_MONTHLY_V4R4 (1993-2016)

  • ECCO_L4_SSH_LLC0090GRID_SNAPSHOT_V4R4 (1993/1/1-2017/1/1, 1st of each month)

  • ECCO_L4_TEMP_SALINITY_LLC0090GRID_SNAPSHOT_V4R4 (1993/1/1-2017/1/1, 1st of each month)

If you haven’t yet been through the download tutorial or used the ecco_download module, it may help you to review that information before downloading the datasets. If you are downloading the snapshots, you may also find the snaps_monthly_textlist function helpful; it generates a file list of only the snapshot files for the 1st of each month, which you can then download using wget.

Prepare environment and load ECCOv4 diagnostic output

Import relevant Python modules

[1]:

import numpy as np
import xarray as xr
import sys
import glob

[2]:

# Suppress warning messages for a cleaner presentation
import warnings
warnings.filterwarnings('ignore')

[3]:

# setting up a dask LocalCluster
from dask.distributed import Client
from dask.distributed import LocalCluster
cluster = LocalCluster()
client = Client(cluster)
client

[3]:

Client

Client-2dad5cfc-ea2e-11ed-9444-c4236006c92f

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:64300/status

Cluster Info

 
[4]:

## Import the ecco_v4_py library into Python
## =========================================
##    If ecco_v4_py is not installed in your local Python library,
##    tell Python where to find it.  The example below adds
##    ecco_v4_py to the user's path if it is stored in the folder
##    ECCOv4-py under the user's home directory

from os.path import join,expanduser,exists,split
user_home_dir = expanduser('~')

sys.path.append(join(user_home_dir,'ECCOv4-py'))

import ecco_v4_py as ecco

[5]:

# Plotting
import matplotlib.pyplot as plt
%matplotlib inline

Add relevant constants

[6]:

# Seawater density (kg/m^3)
rhoconst = 1029
## needed to convert surface mass fluxes to volume fluxes

# Heat capacity (J/kg/K)
c_p = 3994

# Constants for surface heat penetration (from Table 2 of Paulson and Simpson, 1977)
R = 0.62
zeta1 = 0.6
zeta2 = 20.0

Load ecco_grid

[7]:

## Set top-level file directory for the ECCO NetCDF files
## =================================================================

## currently set to ~/Downloads/ECCO_V4r4_PODAAC,
## the default if ecco_podaac_download was used to download dataset granules
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')

[8]:

## Load the model grid
ecco_grid = xr.open_dataset(glob.glob(join(ECCO_dir,'*GEOMETRY*','*.nc'))[0])

Volume

Calculate the volume of each grid cell. This is used when converting advective and diffusive flux convergences and calculating volume-weighted averages.

[9]:

# Volume (m^3)
vol = (ecco_grid.rA*ecco_grid.drF*ecco_grid.hFacC).transpose('tile','k','j','i')

Load monthly snapshots

[10]:

year_start = 1993
year_end = 2017

# open ETAN and THETA snapshots (beginning of each month)
ecco_monthly_SSH = xr.open_mfdataset(join(ECCO_dir,'*SSH*SNAPSHOT*','*-01T*.nc'),\
                                       data_vars='minimal',coords='minimal',compat='override')
ecco_monthly_TS = xr.open_mfdataset(join(ECCO_dir,'*TEMP_SALINITY*SNAPSHOT*','*-01T*.nc'),\
                                       data_vars='minimal',coords='minimal',compat='override')
ecco_monthly_snaps = xr.merge((ecco_monthly_SSH['ETAN'],ecco_monthly_TS['THETA']))

# time mask for snapshots
time_snap_mask = np.logical_and(ecco_monthly_snaps.time.values >= np.datetime64(str(year_start)+'-01-01','ns'),\
                                ecco_monthly_snaps.time.values < np.datetime64(str(year_end)+'-01-02','ns'))

ecco_monthly_snaps = ecco_monthly_snaps.isel(time=time_snap_mask)

[11]:

# 1993-01 (beginning of first month) to 2015-01-01 (end of last month, 2014-12)
print(ecco_monthly_snaps.ETAN.time.isel(time=[0, -1]).values)

['1993-01-01T00:00:00.000000000' '2017-01-01T00:00:00.000000000']

[12]:

# Find the record of the last snapshot
## This is used to defined the exact period for monthly mean data
last_record_date = ecco.extract_yyyy_mm_dd_hh_mm_ss_from_datetime64(ecco_monthly_snaps.time[-1].values)
print(last_record_date)

(2017, 1, 1, 0, 0, 0)

Load monthly mean data

[13]:

## Open ECCO monthly mean variables
ecco_vars_int = xr.open_mfdataset(join(ECCO_dir,'*_OCEAN_3D_TEMPERATURE_FLUX*MONTHLY*','*.nc'),\
                                       data_vars='minimal',coords='minimal',compat='override')
ecco_vars_sfc = xr.open_mfdataset(join(ECCO_dir,'*_HEAT_FLUX*MONTHLY*','*.nc'),\
                                       data_vars='minimal',coords='minimal',compat='override')
ecco_monthly_mean = xr.merge((ecco_vars_int,\
                              ecco_vars_sfc[['TFLUX','oceQsw']]))

# time mask for monthly means
time_mean_mask = np.logical_and(ecco_monthly_mean.time.values >= np.datetime64(str(year_start)+'-01-01','ns'),\
                                ecco_monthly_mean.time.values < np.datetime64(str(year_end)+'-01-01','ns'))

ecco_monthly_mean = ecco_monthly_mean.isel(time=time_mean_mask)

[14]:

# Print first and last time points of the monthly-mean records
print(ecco_monthly_mean.time.isel(time=[0, -1]).values)

['1993-01-16T12:00:00.000000000' '2016-12-16T12:00:00.000000000']

Each monthly mean record is bookended by a snapshot. We should have one more snapshot than monthly mean record.

[15]:

print('Number of monthly mean records: ', len(ecco_monthly_mean.time))
print('Number of monthly snapshot records: ', len(ecco_monthly_snaps.time))

Number of monthly mean records:  288
Number of monthly snapshot records:  289

[16]:

# Drop superfluous coordinates (We already have them in ecco_grid)
ecco_monthly_mean = ecco_monthly_mean.reset_coords(drop=True)

Merge dataset of monthly mean and snapshots data

Merge the two datasets to put everything into one single dataset

[17]:

ds = xr.merge([ecco_monthly_mean,
               ecco_monthly_snaps.rename({'time':'time_snp','ETAN':'ETAN_snp', 'THETA':'THETA_snp'})])

Create the xgcm ‘grid’ object

The xgcm ‘grid’ object is used to calculate the flux divergences across different tiles of the lat-lon-cap grid and the time derivatives from THETA snapshots

[18]:

# Change time axis of the snapshot variables
ds.time_snp.attrs['c_grid_axis_shift'] = 0.5

[19]:

grid = ecco.get_llc_grid(ds)

Number of seconds in each month

The xgcm grid object includes information on the time axis, such that we can use it to get \Delta t, which is the time span between the beginning and end of each month (in seconds).

[20]:

delta_t = grid.diff(ds.time_snp, 'T', boundary='fill', fill_value=np.nan)

# Convert to seconds
delta_t = delta_t.astype('f4') / 1e9

Calculate total tendency of \theta (G^{\theta}_\textrm{total})

We calculate the monthly-averaged time tendency of THETA by differencing monthly THETA snapshots. Remember that we need to include a scaling factor due to the nonlinear free surface formulation. Thus, we need to use snapshots of both ETAN and THETA to evaluate s^*\theta.

[21]:

# Calculate the s*theta term
sTHETA = ds.THETA_snp*(1+ds.ETAN_snp/ecco_grid.Depth)

[22]:

# Total tendency (psu/s)
G_total = sTHETA.diff(dim='time_snp')/np.expand_dims(delta_t.values,axis=(1,2,3,4))

# re-assign and rename time coordinate
G_total = G_total.rename({'time_snp':'time'})
G_total = G_total.assign_coords({'time':delta_t.time.values})

[23]:

G_total

[23]:

<xarray.DataArray (time: 288, k: 50, tile: 13, j: 90, i: 90)>
dask.array<truediv, shape=(288, 50, 13, 90, 90), dtype=float32, chunksize=(1, 50, 13, 90, 90), chunktype=numpy.ndarray>
Coordinates:
  * i        (i) int32 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
  * j        (j) int32 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
  * k        (k) int32 0 1 2 3 4 5 6 7 8 9 10 ... 40 41 42 43 44 45 46 47 48 49
  * tile     (tile) int32 0 1 2 3 4 5 6 7 8 9 10 11 12
  * time     (time) datetime64[ns] 1993-01-16T12:00:00 ... 2016-12-16T12:00:00
    XC       (tile, j, i) float32 -111.6 -111.3 -110.9 ... -99.42 -105.6 -111.9
    YC       (tile, j, i) float32 -88.24 -88.38 -88.52 ... -88.03 -88.08 -88.1
    Z        (k) float32 dask.array<chunksize=(50,), meta=np.ndarray>

Note: Unlike the monthly snapshots ETAN_snp and THETA_snp, the resulting data array G_total has now the same time values as the time-mean fields (middle of the month).

Plot the time-mean \partial \theta / \partial t, total \Delta \theta, and one example \partial \theta / \partial t field

Time-mean \partial \theta / \partial t

The time-mean \partial \theta / \partial t (i.e., \overline{G^{\theta}_\textrm{total}}), is given by

\overline{G^{\theta}_\textrm{total}} = \sum_{i=1}^{nm} w_i G^{\theta}_\textrm{total}

with \sum_{i=1}^{nm} w_i = 1 and nm=number of months

[24]:

# The weights are just the number of seconds per month divided by total seconds
month_length_weights = delta_t / delta_t.sum()

[25]:

# The weighted mean weights by the length of each month (in seconds)
G_total_mean = (G_total*month_length_weights).sum('time').compute()

[26]:

plt.figure(figsize=(15,15))

for idx, k in enumerate([0,10,25]):
    p = ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_total_mean.isel(k=k),show_colorbar=True,
                                      cmap='RdBu_r', user_lon_0=-67, dx=2, dy=2, subplot_grid=[3,1,idx+1]);
    p[1].set_title(r'$\overline{G^\theta_{total}}$ at z = %i m (k = %i) [$^\circ$C s$^{-1}$]'\
                   %(np.round(-ecco_grid.Z[k].values),k), fontsize=16)

-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_47_1.png

Total \Delta \theta

How much did THETA change over the analysis period?

[27]:

# The number of seconds in the entire period
seconds_in_entire_period = \
    float(ds.time_snp[-1] - ds.time_snp[0])/1e9
print ('seconds in analysis period: ', seconds_in_entire_period)

# which is also the sum of the number of seconds in each month
print('Sum of seconds in each month ', delta_t.sum().values)

seconds in analysis period:  757382400.0
Sum of seconds in each month  757382400.0

[28]:

THETA_delta = G_total_mean*seconds_in_entire_period

[29]:

plt.figure(figsize=(15,5));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              THETA_delta.isel(k=0),show_colorbar=True,\
                              cmin=-4, cmax=4, \
                              cmap='RdBu_r', user_lon_0=-67, dx=0.2, dy=0.2);
plt.title(r'Predicted $\Delta \theta$ at the sea surface [$^\circ$C] from $\overline{G^\theta_{total}}$',fontsize=16);

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_51_1.png

We can sanity check the total THETA change that we found by multipling the time-mean THETA tendency with the number of seconds in the simulation by comparing that with the difference in THETA between the end of the last month and start of the first month.

[30]:

THETA_delta_method_2 = ds.THETA_snp.isel(time_snp=-1) - ds.THETA_snp.isel(time_snp=0)

[31]:

plt.figure(figsize=(15,5));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              THETA_delta_method_2.isel(k=0),show_colorbar=True,\
                              cmin=-4, cmax=4, \
                              cmap='RdBu_r', user_lon_0=-67, dx=0.2, dy=0.2);
plt.title(r'Actual $\Delta \theta$ [$^\circ$C]', fontsize=16);

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_54_1.png

Example G^\theta_{total} field at a particular time

[32]:

# get an array of YYYY, MM, DD, HH, MM, SS for
#dETAN_dT_perSec at time index 100
curr_t_ind = 100
tmp = str(G_total.time.values[curr_t_ind])
print(tmp)

2001-05-16T12:00:00.000000000

[33]:

plt.figure(figsize=(15,5));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_total.isel(time=curr_t_ind,k=0), show_colorbar=True,
                              cmap='RdBu_r', user_lon_0=-67, dx=0.2, dy=0.2);

plt.title(r'$G^\theta_{total}$ at the sea surface [$^\circ$C s$^{-1}$] during ' +
          str(tmp)[0:4] +'/' + str(tmp)[5:7], fontsize=16);

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_57_1.png

For any given month the time rate of change of THETA is strongly dependent on the season. In the above we are looking at May 2001. We see positive THETA tendency in the northern hemisphere and cooling in the southern hemisphere.

Calculate tendency due to advective convergence (G^{\theta}_\textrm{advection})

Horizontal convergence of advective heat flux

The relevant fields from the diagnostic output here are - ADVx_TH: U Component Advective Flux of Potential Temperature (degC m^3/s) - ADVy_TH: V Component Advective Flux of Potential Temperature (degC m^3/s)

The xgcm grid object is then used to take the convergence of the horizontal heat advection.

Note: when using at least one recent version of xgcm (v0.8.1), errors were triggered when calling diff_2d_vector below. If you have this issue, as a temporary fix you can locate the xgcm directory

import xgcm
print(xgcm.__path__)

and then substitute padding.py and grid_ufunc.py with the files linked below:

padding.py

grid_ufunc.py

Alternatively, you can revert xgcm to an earlier version (e.g., v0.6.0) that should not have this issue.

[34]:

# Set fluxes on land to zero (instead of NaN)
ds['ADVx_TH'] = ds.ADVx_TH.where(ecco_grid.hFacW.values > 0,0)
ds['ADVy_TH'] = ds.ADVy_TH.where(ecco_grid.hFacS.values > 0,0)

# compute horizontal components of divergence
ADVxy_diff = grid.diff_2d_vector({'X' : ds.ADVx_TH, 'Y' : ds.ADVy_TH}, boundary = 'fill')

# Convergence of horizontal advection (degC m^3/s)
adv_hConvH = (-(ADVxy_diff['X'] + ADVxy_diff['Y']))

Vertical convergence of advective heat flux

The relevant field from the diagnostic output is - ADVr_TH: Vertical Advective Flux of Potential Temperature (degC m^3/s)

[35]:

# Set fluxes on land to zero (instead of NaN)
ds['ADVr_TH'] = ds.ADVr_TH.where(ecco_grid.hFacC.values > 0,0)

# Load monthly averages of vertical advective flux
ADVr_TH = ds.ADVr_TH.transpose('time','tile','k_l','j','i')

Note: For ADVr_TH, DFrE_TH and DFrI_TH, we need to make sure that sequence of dimensions are consistent. When loading the fields use .transpose('time','tile','k_l','j','i'). Otherwise, the divergences will be not correct (at least for tile = 12).

[36]:

# Convergence of vertical advection (degC m^3/s)
adv_vConvH = grid.diff(ADVr_TH, 'Z', boundary='fill')

Note: In case of the volume budget (and salinity conservation), the surface forcing (oceFWflx) is already included at the top level (k_l = 0) in WVELMASS. Thus, to keep the surface forcing term explicitly represented, one needs to zero out the values of WVELMASS at the surface so as to avoid double counting (see the Volume budget closure tutorial). This is not the case for the heat budget. ADVr_TH does not include the sea surface forcing. Thus, the vertical advective flux (at the air-sea interface) should not be zeroed out.

Total convergence of advective flux (G^{\theta}_\textrm{advection})

We can get the total convergence by simply adding the horizontal and vertical component.

[37]:

# Sum horizontal and vertical convergences and divide by volume (degC/s)
G_advection = (adv_hConvH + adv_vConvH)/vol

Plot the time-mean G^{\theta}_\textrm{advection}

[38]:

G_advection_mean = (G_advection*month_length_weights).sum('time').compute()

[39]:

plt.figure(figsize=(15,15))

for idx, k in enumerate([0,1,25]):
    p = ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_advection_mean.isel(k=k),show_colorbar=True,
                                      cmin=-1e-6, cmax=1e-6, cmap='RdBu_r', user_lon_0=-67, dx=2, dy=2,
                                      subplot_grid=[3,1,idx+1]);
    p[1].set_title(r'$\overline{G^\theta_{advection}}$ at z = %i m (k = %i) [$^\circ$C s$^{-1}$]'\
                   %(np.round(-ecco_grid.Z[k].values),k), fontsize=16)

-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_70_1.png

Example G^{\theta}_\textrm{advection} field at a particular time

[40]:

curr_t_ind = 100
tmp = str(G_advection.time.values[curr_t_ind])
print(tmp)

2001-05-16T12:00:00.000000000

[41]:

plt.figure(figsize=(15,5));

ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_advection.isel(time=curr_t_ind,k=k), show_colorbar=True,
                              cmin=-1e-6, cmax=1e-6, cmap='RdBu_r', user_lon_0=-67, dx=0.2, dy=0.2)
plt.title(r'$G^\theta_{advection}$ at the sea surface [$^\circ$C s$^{-1}$] during ' +
          str(tmp)[0:4] +'/' + str(tmp)[5:7], fontsize=16)
plt.show()

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_73_1.png

Calculate tendency due to diffusive convergence (G^{\theta}_\textrm{diffusion})

Horizontal convergence of advective heat flux

The relevant fields from the diagnostic output here are - DFxE_TH: U Component Diffusive Flux of Potential Temperature (degC m^3/s) - DFyE_TH: V Component Diffusive Flux of Potential Temperature (degC m^3/s)

As with advective fluxes, we use the xgcm grid object to calculate the convergence of horizontal heat diffusion.

[42]:

# Set fluxes on land to zero (instead of NaN)
ds['DFxE_TH'] = ds.DFxE_TH.where(ecco_grid.hFacW.values > 0,0)
ds['DFyE_TH'] = ds.DFyE_TH.where(ecco_grid.hFacS.values > 0,0)

DFxyE_diff = grid.diff_2d_vector({'X' : ds.DFxE_TH, 'Y' : ds.DFyE_TH}, boundary = 'fill')

# Convergence of horizontal diffusion (degC m^3/s)
dif_hConvH = (-(DFxyE_diff['X'] + DFxyE_diff['Y']))

Vertical convergence of advective heat flux

The relevant fields from the diagnostic output are - DFrE_TH: Vertical Diffusive Flux of Potential Temperature (Explicit part) (degC m^3/s) - DFrI_TH: Vertical Diffusive Flux of Potential Temperature (Implicit part) (degC m^3/s) > Note: Vertical diffusion has both an explicit (DFrE_TH) and an implicit (DFrI_TH) part.

[43]:

# Set land values to zero (not NaN)
ds['ADVr_TH'] = ds.ADVr_TH.where(ecco_grid.hFacC.values > 0,0)

# Set fluxes on land to zero (instead of NaN)
ds['DFrE_TH'] = ds.DFrE_TH.where(ecco_grid.hFacC.values > 0,0)
ds['DFrI_TH'] = ds.DFrI_TH.where(ecco_grid.hFacC.values > 0,0)

# Load monthly averages of vertical diffusive fluxes
DFrE_TH = ds.DFrE_TH.transpose('time','tile','k_l','j','i')
DFrI_TH = ds.DFrI_TH.transpose('time','tile','k_l','j','i')

# Convergence of vertical diffusion (degC m^3/s)
dif_vConvH = grid.diff(DFrE_TH, 'Z', boundary='fill') + grid.diff(DFrI_TH, 'Z', boundary='fill')

Total convergence of diffusive flux (G^{\theta}_\textrm{diffusion})

[44]:

# Sum horizontal and vertical convergences and divide by volume (degC/s)
G_diffusion = (dif_hConvH + dif_vConvH)/vol

Plot the time-mean G^{\theta}_\textrm{diffusion}

[45]:

G_diffusion_mean = (G_diffusion*month_length_weights).sum('time').compute()

[46]:

plt.figure(figsize=(15,15))

for idx, k in enumerate([0,1,25]):
    p = ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_diffusion_mean.isel(k=k),show_colorbar=True,
                                      cmin=-3e-6, cmax=3e-6, cmap='RdBu_r', user_lon_0=-67, dx=2, dy=2,
                                      subplot_grid=[3,1,idx+1]);
    p[1].set_title(r'$\overline{G^\theta_{diffusion}}$ at z = %i m (k = %i) [$^\circ$C s$^{-1}$]'\
                   %(np.round(-ecco_grid.Z[k].values),k), fontsize=16)

-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_82_1.png

Example G^{\theta}_\textrm{diffusion} field at a particular time

[47]:

curr_t_ind = 100
tmp = str(G_diffusion.time.values[curr_t_ind])
print(tmp)

2001-05-16T12:00:00.000000000

[48]:

plt.figure(figsize=(15,5));

ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_diffusion.isel(time=curr_t_ind,k=0),show_colorbar=True,
                              cmin=-3e-6, cmax=3e-6, cmap='RdBu_r', user_lon_0=-67, dx=0.2, dy=0.2)
plt.title(r'$G^\theta_{diffusion}$ at the sea surface [$^\circ$C s$^{-1}$] during ' +
          str(tmp)[0:4] +'/' + str(tmp)[5:7], fontsize=16)
plt.show()

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_85_1.png

Calculate tendency due to forcing (G^{\theta}_\textrm{forcing})

Finally, we evaluate the local forcing term due to surface heat and geothermal fluxes.

Surface heat flux

For the surface contribution, there are two relevant model diagnostics: - TFLUX: total heat flux (match heat-content variations) (W/m^2) - oceQsw: net Short-Wave radiation (+=down) (W/m^2)

Defining terms needed for evaluating surface heat forcing

[49]:

Z = ecco_grid.Z.compute()
RF = np.concatenate([ecco_grid.Zp1.values[:-1],[np.nan]])

Note: Z and Zp1 are used in deriving surface heat penetration. MATLAB code uses RF from mygrid structure.

[50]:

q1 = R*np.exp(1.0/zeta1*RF[:-1]) + (1.0-R)*np.exp(1.0/zeta2*RF[:-1])
q2 = R*np.exp(1.0/zeta1*RF[1:]) + (1.0-R)*np.exp(1.0/zeta2*RF[1:])

[51]:

# Correction for the 200m cutoff
zCut = np.where(Z < -200)[0][0]
q1[zCut:] = 0
q2[zCut-1:] = 0

[52]:

# Create xarray data arrays
q1 = xr.DataArray(q1,coords=[Z.k],dims=['k'])
q2 = xr.DataArray(q2,coords=[Z.k],dims=['k'])

Compute vertically penetrating flux

Given the penetrating nature of the shortwave term, to properly evaluate the local forcing term, oceQsw must be removed from TFLUX (which contains the net latent, sensible, longwave, and shortwave contributions) and redistributed vertically.

[53]:

## Land masks
# Make copy of hFacC
mskC = ecco_grid.hFacC.copy(deep=True).compute()

# Change all fractions (ocean) to 1. land = 0
mskC.values[mskC.values>0] = 1

[54]:

# Shortwave flux below the surface (W/m^2)
forcH_subsurf = ((q1*(mskC==1)-q2*(mskC.shift(k=-1)==1))*ds.oceQsw).transpose('time','tile','k','j','i')

[55]:

# Surface heat flux (W/m^2)
forcH_surf = ((ds.TFLUX - (1-(q1[0]-q2[0]))*ds.oceQsw)\
              *mskC[0]).transpose('time','tile','j','i').assign_coords(k=0).expand_dims('k')

[56]:

# Full-depth sea surface forcing (W/m^2)
forcH = xr.concat([forcH_surf,forcH_subsurf[:,:,1:]], dim='k').transpose('time','tile','k','j','i')

Geothermal flux

The geothermal flux contribution is not accounted for in any of the standard model diagnostics provided as output. Rather, this term, which is time invariant, is provided in the input file geothermalFlux.bin contained in the ancillary data archive. This archive can be downloaded from PO.DAAC, though accessing it there requires downloading a much larger tarball of files (~192 GB). So, for the time being, the geothermalFlux.bin file is also stored on the tutorial Github and can be downloaded here.

Note: The code cell below assumes geothermalFlux.bin has been placed in ~/Downloads after downloading. Change the directory as needed.

[57]:

# Load the geothermal heat flux using the routine 'read_llc_to_tiles'.
geoflx_dir = join(user_home_dir,'Downloads')
geoflx = ecco.read_llc_to_tiles(geoflx_dir, 'geothermalFlux.bin')

load_binary_array: loading file C:\Users\adelman\Downloads\geothermalFlux.bin
load_binary_array: data array shape  (1170, 90)
load_binary_array: data array type  >f4
llc_compact_to_faces: dims, llc  (1170, 90) 90
llc_compact_to_faces: data_compact array type  >f4
llc_faces_to_tiles: data_tiles shape  (13, 90, 90)
llc_faces_to_tiles: data_tiles dtype  >f4

The geothermal flux dataset needs to be saved as an xarray data array with the same format as the model output.

[58]:

# Convert numpy array to an xarray DataArray with matching dimensions as the monthly mean fields
geoflx_llc = xr.DataArray(geoflx,coords={'tile': ecco_monthly_mean.tile.values,
                                         'j': ecco_monthly_mean.j.values,
                                         'i': ecco_monthly_mean.i.values},dims=['tile','j','i'])

[59]:

plt.figure(figsize=(15,5));

ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, geoflx_llc,show_colorbar=True,cmap='magma',
                              user_lon_0=-67, dx=0.2, dy=0.2)
plt.title(r'Geothermal heat flux [W m$^{-2}$]', fontsize=16)
plt.show()

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_104_1.png

Geothermal flux needs to be a three dimensional field since the sources are distributed along the ocean floor at various depths. This requires a three dimensional mask.

[60]:

# Create 3d bathymetry mask
mskC_shifted = mskC.shift(k=-1)

mskC_shifted.values[-1,:,:,:] = 0
mskb = mskC - mskC_shifted

# Create 3d field of geothermal heat flux
geoflx3d = geoflx_llc * mskb.transpose('k','tile','j','i')
GEOFLX = geoflx3d.transpose('k','tile','j','i')
GEOFLX.attrs = {'standard_name': 'GEOFLX','long_name': 'Geothermal heat flux','units': 'W/m^2'}

Total forcing (G^{\theta}_\textrm{forcing})

[61]:

# Add geothermal heat flux to forcing field and convert from W/m^2 to degC/s
G_forcing = ((forcH + GEOFLX)/(rhoconst*c_p))/(ecco_grid.hFacC*ecco_grid.drF)

Plot the time-mean G^{\theta}_\textrm{forcing}

[62]:

G_forcing_mean = (G_forcing*month_length_weights).sum('time')

[63]:

plt.figure(figsize=(15,15))

for idx, k in enumerate([0,1,25]):
    p = ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_forcing_mean.isel(k=k),show_colorbar=True,
                                      cmin=-3e-6, cmax=3e-6, cmap='RdBu_r', user_lon_0=-67, dx=2, dy=2,
                                      subplot_grid=[3,1,idx+1]);
    p[1].set_title(r'$\overline{G^\theta_{forcing}}$ at z = %i m (k = %i) [$^\circ$C s$^{-1}$]'\
                   %(np.round(-ecco_grid.Z[k].values),k), fontsize=16)

-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
-178.99657534246575 111.99657534246575
-180.0 113.0
-89.0 89.0
-90.0 90.0
114.01516136363637 178.98484863636364
113.00001 180.0
-89.0 89.0
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_111_1.png

\overline{G^\theta_{forcing}} is focused at the sea surface and much smaller (essentially zero) at depth. \overline{G^\theta_{forcing}} is negative for most of the ocean (away from the equator). The spatial pattern in the surface forcing is the same as for diffusion but with opposite sign (see maps for \overline{G^\theta_{diffusion}} above). This makes sense as forcing is to a large extent balanced by diffusion within the mixed layer.

Example G^{\theta}_\textrm{forcing} field at a particular time

[64]:

curr_t_ind = 100
tmp = str(G_advection.time.values[curr_t_ind])
print(tmp)

2001-05-16T12:00:00.000000000

[65]:

plt.figure(figsize=(15,5));

ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_forcing.isel(time=curr_t_ind,k=0),show_colorbar=True,
                              cmin=-5e-6, cmax=5e-6, cmap='RdBu_r', user_lon_0=-67, dx=0.2, dy=0.2)
plt.title(r'$G^\theta_{forcing}$ at the sea surface [$^\circ$C s$^{-1}$] during ' +
          str(tmp)[0:4] +'/' + str(tmp)[5:7], fontsize=16)
plt.show()

-179.9 112.9
-180.0 113.0
-89.9 89.9
-90.0 90.0
113.10030938622755 179.89970061377244
113.00001 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_115_1.png

Save to dataset

Now that we have all the terms evaluated, let’s save them to a dataset. Here are two examples: - Zarr is a new format that is used for cloud storage. - Netcdf is the more traditional format that most people are familiar with.

When saving this heat budget dataset, the zarr file is ~15 GB, while the NetCDF file is ~53 GB. So zarr can be more efficient for storage.

Add all variables to a new dataset

[66]:

varnames = ['G_total','G_advection','G_diffusion','G_forcing']

ds = xr.Dataset(data_vars={})
for varname in varnames:
    ds[varname] = globals()[varname].chunk(chunks={'time':1,'tile':13,'k':50,'j':90,'i':90})

[67]:

# Add surface forcing (degC/s)
ds['Qnet'] = ((forcH /(rhoconst*c_p))\
              /(ecco_grid.hFacC*ecco_grid.drF)).chunk(chunks={'time':1,'tile':13,'k':50,'j':90,'i':90})

[68]:

# Add shortwave penetrative flux (degC/s)
#Since we only are interested in the subsurface heat flux we need to zero out the top cell
SWpen = ((forcH_subsurf /(rhoconst*c_p))/(ecco_grid.hFacC*ecco_grid.drF)).where(forcH_subsurf.k>0).fillna(0.)
ds['SWpen'] = SWpen.where(ecco_grid.hFacC>0).chunk(chunks={'time':1,'tile':13,'k':50,'j':90,'i':90})

Note: Qnet and SWpen are included in G_forcing and are not necessary to close the heat budget.

[69]:

ds.time.encoding = {}
ds = ds.reset_coords(drop=True)

Save to zarr

[70]:

from dask.diagnostics import ProgressBar

[71]:

# save_dir is set to ~/Downloads below;
# change if you want to save somewhere else
save_dir = join(user_home_dir,'Downloads')
with ProgressBar():
    ds.to_zarr(join(save_dir,'eccov4r4_budg_heat'))

Save to netcdf

[72]:

with ProgressBar():
    ds.to_netcdf(join(save_dir,'eccov4r4_budg_heat.nc'), format='NETCDF4')

Load budget variables from file

After having saved the budget terms to file, we can load the dataset like this

[73]:

# Load terms from zarr dataset
G_total = xr.open_zarr(join(save_dir,'eccov4r4_budg_heat')).G_total
G_advection = xr.open_zarr(join(save_dir,'eccov4r4_budg_heat')).G_advection
G_diffusion = xr.open_zarr(join(save_dir,'eccov4r4_budg_heat')).G_diffusion
G_forcing = xr.open_zarr(join(save_dir,'eccov4r4_budg_heat')).G_forcing
Qnet = xr.open_zarr(join(save_dir,'eccov4r4_budg_heat')).Qnet
SWpen = xr.open_zarr(join(save_dir,'eccov4r4_budg_heat')).SWpen

Or if you saved it as a netcdf file:

# Load terms from netcdf file
G_budget = xr.open_mfdataset(join(save_dir,'eccov4r4_budg_heat.nc'))
G_total_tendency = G_budget.G_total_tendency
G_advection = G_budget.G_advection
G_diffusion = G_budget.G_diffusion
G_forcing = G_budget.G_forcing
Qnet = G_budget.Qnet

Comparison between LHS and RHS of the budget equation

[74]:

# Total convergence
ConvH = G_advection + G_diffusion

[75]:

# Sum of terms in RHS of equation
rhs = ConvH + G_forcing

Map of residuals

[76]:

res = ((rhs-G_total).sum(dim='k')*month_length_weights).sum(dim='time').compute()

[77]:

plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, res,
                              cmin=-5e-12, cmax=5e-12, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Residual $\partial \theta / \partial t$ [$^\circ$C s$^{-1}$]: RHS - LHS', fontsize=16)
plt.show()

-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_135_1.png

The residual (summed over depth and time) is essentially zero everywhere. What if we omit the geothermal heat flux?

[78]:

# Residual when omitting geothermal heat flux
res_geo = ((ConvH + Qnet - G_total).sum(dim='k')*month_length_weights).sum(dim='time').compute()

[79]:

plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, res_geo,
                              cmin=-1e-9, cmax=1e-9, show_colorbar=True, cmap='RdBu_r', dx=0.2, dy=0.2)
plt.title(r'Residual due to omitting geothermal heat [$^\circ$C s$^{-1}$] ', fontsize=16)
plt.show()

-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_138_1.png

We see that the contribution from geothermal flux in the heat budget is well above the residual (by three orders of magnitude).

[80]:

# Residual when omitting shortwave penetrative heat flux
res_sw = (rhs-SWpen-G_total).sum(dim='k').sum(dim='time').compute()

[81]:

plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, res_sw,
                              cmin=-5e-4, cmax=5e-4, show_colorbar=True, cmap='RdBu_r', dx=0.2, dy=0.2)
plt.title(r'Residual due to omitting shortwave penetrative heat flux [$^\circ$C s$^{-1}$] ', fontsize=16)
plt.show()

-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
_images/ECCO_v4_Heat_budget_closure_141_1.png

In terms of subsurface heat fluxes, shortwave penetration represents a much larger heat flux compared to geothermal heat flux (by around three orders of magnitude).

Histogram of residuals

We can look at the distribution of residuals to get a little more confidence.

[82]:

tmp = np.abs(res).values.ravel()

[83]:

plt.figure(figsize=(10,3));

plt.hist(tmp[tmp > 0],np.linspace(0, 2.e-12, 200));
plt.grid()
_images/ECCO_v4_Heat_budget_closure_145_0.png

Summing residuals vertically and temporally yields < 10^{-12} ^\circC s^{-1} for most grid points.

Heat budget closure through time

Global average budget closure

Another way of demonstrating heat budget closure is to show the global spatially-averaged THETA tendency terms

[84]:

# Volume (m^3)
vol = (ecco_grid.rA*ecco_grid.drF*ecco_grid.hFacC).transpose('tile','k','j','i')

# Take volume-weighted mean of these terms
tmp_a=(G_total*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_b=(G_advection*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_c=(G_diffusion*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_d=(G_forcing*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_e=(rhs*vol).sum(dim=('k','i','j','tile'))/vol.sum()

# Result is five time series
tmp_a.dims

[84]:

('time',)

[85]:

fig, axs = plt.subplots(2, 2, figsize=(14,8))

plt.sca(axs[0,0])
tmp_a.plot(color='k',lw=2)
tmp_e.plot(color='grey')
axs[0,0].set_title(r'a. $G^\theta_{total}$ (black) / RHS (grey) [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()

plt.sca(axs[0,1])
tmp_b.plot(color='r')
axs[0,1].set_title(r'b. $G^\theta_{advection}$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()

plt.sca(axs[1,0])
tmp_c.plot(color='orange')
axs[1,0].set_title(r'c. $G^\theta_{diffusion}$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()

plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title(r'd. $G^\theta_{forcing}$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Global Heat Budget', fontsize=16);
_images/ECCO_v4_Heat_budget_closure_149_0.png

When averaged over the entire ocean the ocean heat transport terms (G^\theta_\textrm{advection} and G^\theta_\textrm{diffusion}) have no net impact on G^\theta_\textrm{total} (i.e., \partial \theta / \partial t). This makes sense because G^\theta_\textrm{advection} and G^\theta_\textrm{diffusion} can only redistributes heat. Globally, \theta can only change via G^\theta_\textrm{forcing}.

Local heat budget closure

Locally we expect that heat divergence can impact \theta. This is demonstrated for a single grid point.

[86]:

# Pick any set of indices (tile, k, j, i) corresponding to an ocean grid point
t,k,j,i = (6,10,40,29)
print(t,k,j,i)

6 10 40 29

[87]:

tmp_a = G_total.isel(tile=t,k=k,j=j,i=i)
tmp_b = G_advection.isel(tile=t,k=k,j=j,i=i)
tmp_c = G_diffusion.isel(tile=t,k=k,j=j,i=i)
tmp_d = G_forcing.isel(tile=t,k=k,j=j,i=i)
tmp_e = rhs.isel(tile=t,k=k,j=j,i=i)

fig, axs = plt.subplots(2, 2, figsize=(14,8))

plt.sca(axs[0,0])
tmp_a.plot(color='k',lw=2)
tmp_e.plot(color='grey')
axs[0,0].set_title(r'a. $G^\theta_{total}$ (black) / RHS (grey) [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()

plt.sca(axs[0,1])
tmp_b.plot(color='r')
axs[0,1].set_title(r'b. $G^\theta_{advection}$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()

plt.sca(axs[1,0])
tmp_c.plot(color='orange')
axs[1,0].set_title(r'c. $G^\theta_{diffusion}$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()

plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title(r'd. $G^\theta_{forcing}$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Heat Budget for one grid point (tile = %i, k = %i, j = %i, i = %i)'%(t,k,j,i), fontsize=16);
_images/ECCO_v4_Heat_budget_closure_152_0.png

Indeed, the heat divergence terms do contribute to \theta variations at a single point. Local heat budget closure is also confirmed at this grid point as we see that the sum of terms on the RHS (grey line) equals the LHS (black line).

For the Arctic grid point, there is a clear seasonal cycles in both G^\theta_\textrm{total}, G^\theta_\textrm{diffusion} and G^\theta_\textrm{forcing}. The seasonal cycle in G^\theta_\textrm{forcing} seems to be the reverse of G^\theta_\textrm{total} and G^\theta_\textrm{diffusion}.

[88]:

plt.figure(figsize=(10,6));
tmp_a.groupby('time.month').mean('time').plot(color='k',lw=3)
tmp_b.groupby('time.month').mean('time').plot(color='r')
tmp_c.groupby('time.month').mean('time').plot(color='orange')
tmp_d.groupby('time.month').mean('time').plot(color='b')
tmp_e.groupby('time.month').mean('time').plot(color='grey')
plt.ylabel(r'$\partial\theta$/$\partial t$ [$^\circ$C s$^{-1}$]', fontsize=12)
plt.grid()
plt.title('Climatological seasonal cycles', fontsize=14)
plt.show()
_images/ECCO_v4_Heat_budget_closure_154_0.png

The mean seasonal cycle of the total is a balance between advection and diffusion. However, this is likely depth-dependent. How does the balance look across the upper 200 meter at that location?

Time-mean vertical profiles

[89]:

fig = plt.subplots(1, 2, sharey=True, figsize=(12,7))

plt.subplot(1, 2, 1)
plt.plot(G_total.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
         lw=4, color='black', marker='.', label=r'$G^\theta_{total}$ (LHS)')

plt.plot(G_advection.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
         lw=2, color='red', marker='.', label=r'$G^\theta_{advection}$')

plt.plot(G_diffusion.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
         lw=2, color='orange', marker='.', label=r'$G^\theta_{diffusion}$')

plt.plot(G_forcing.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
         lw=2, color='blue', marker='.', label=r'$G^\theta_{forcing}$')
plt.plot(rhs.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z, lw=1, color='grey', marker='.', label='RHS')
plt.xlabel(r'$\partial\theta$/$\partial t$ [$^\circ$C s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.ylabel('Depth (m)', fontsize=14)
plt.gca().tick_params(axis='both', which='major', labelsize=12)
plt.legend(loc='lower left', frameon=False, fontsize=12)

plt.subplot(1, 2, 2)
plt.plot(G_total.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
         lw=4, color='black', marker='.', label=r'$G^\theta_{total}$ (LHS)')
plt.plot(rhs.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z, lw=1, color='grey', marker='.', label='RHS')
plt.setp(plt.gca(), 'yticklabels',[])
plt.xlabel(r'$\partial\theta$/$\partial t$ [$^\circ$C s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.show()
_images/ECCO_v4_Heat_budget_closure_157_0.png

Balance between surface forcing and diffusion in the top layers. Balance between advection and diffusion at depth.