Case: b.e30_alpha08b.B1850C_LTso.ne30_t232_wgx3.sfc_test01

caseroot: /glade/work/chengz/cesm.cases/sfc_fluxes_check/b.e30_alpha08b.B1850C_LTso.ne30_t232_wgx3.sfc_test01

Purposes

Surface fluxes check for the B case.

MOM6 diagnostic boundary fluxes of scalars and their global budgets

This notebook is meant to work with MOM6 run via the CESM framework

b.e30_alpha08b.B1850C_LTso.ne30_t232_wgx3.sfc_test01.

Results from this notebook:

1. Maps of surface boundary fluxes of water, heat, and salt crossing into the liquid seawater in MOM6;

2. Computation of self-consistency checks, including global heat, salt, and mass budgets to verify that the model is conserving scalars over the global/regional domain.

Hopes for the use of this notebook:

1. To provide a starting point to document boundary fluxes of scalar fields;

2. To teach MOM6 users about the boundary fluxes, their patterns, units, and sign conventions;

3. To perform self-consistency checks to ensure the model is conserving scalar fields;

4. To illustrate a self-contained notebook of use for MOM6 analysis.

This notebook was originally developed at NOAA/GFDL, and it is provided freely to the MOM6 community. GFDL scientists developing MOM6 make extensive use of Python for diagnostics. We solicit modifications/fixes that are useful to the MOM6 community.

import matplotlib.pyplot as plt
import xarray as xr
import numpy as np

path    = '/glade/u/home/chengz/scratch/b.e30_alpha08b.B1850C_LTso.ne30_t232_wgx3.sfc_test01/run/'
static  = xr.open_dataset(path+'ocean.mom6.static.nc')
forcing = xr.open_dataset(path+'ocean.mom6.frc.0001.nc', decode_times=False)
surface = xr.open_dataset(path+'ocean.mom6.sfc.0001.nc', decode_times=False)

cp   = 3992.0
rho0 = 1035.0
n    = 3
print('Variables saved:')
for v in forcing.data_vars: print(v)
for v in surface.data_vars: print(v)
for v in static.data_vars: print(v)

Variables saved:
net_massin
net_massout
evap
seaice_melt
lrunoff
frunoff
lrunoff_glc
frunoff_glc
lprec
fprec
vprec
PRCmE
internal_heat
net_heat_surface
net_heat_coupler
LwLatSens
SW
frazil
seaice_melt_heat
Heat_PmE
heat_content_surfwater
heat_content_massin
heat_content_massout
heat_content_evap
heat_content_lprec
heat_content_fprec
heat_content_vprec
heat_content_lrunoff
heat_content_frunoff
heat_content_frunoff_glc
heat_content_lrunoff_glc
heat_content_cond
salt_flux
salt_flux_in
average_T1
average_T2
average_DT
time_bounds
SST
mass_wt
temp_int
salt_int
geolon
geolat
area_t
wet

# This section fills the fields used in this notebook.

#--------------------------------------------------------------
# geometric factors 
lon  = static['geolon']
lat  = static['geolat']
wet  = static['wet']
area = static['area_t'] * wet


#--------------------------------------------------------------
# time in seconds (time coordinate values are in days)
time = surface['time'].values * 86400.


#--------------------------------------------------------------
# sea surface temperature 
sst = surface['SST'].isel(time=n)


#--------------------------------------------------------------
# \int \rho \, dz \, (1,Theta,S), with \rho = \rho_{0} for Boussinesq.  
mass_wt = surface['mass_wt']
tomint  = surface['temp_int']
somint  = surface['salt_int']

#--------------------------------------------------------------

# This case has a daily output frequency.
# The variables above are snapshots at the end of the day (time stamp at 12 am)
# The variables below are daily means over the previous day (time stamp at 12 pm)

#--------------------------------------------------------------
# mass flux of water crossing ocean surface [kg/(m^2 s)]
# positive values indicate mass entering ocean; 
# negative values indicate mass leaving ocean. 

# net mass flux entering ocean 
net_massin = forcing['net_massin'].isel(time=n)

# net mass flux leaving ocean
net_massout = forcing['net_massout'].isel(time=n)

# evaporation (negative) and condensation (positive)
evap = forcing['evap'].isel(time=n)

# meltw, fresh water from melting/forming sea ice
meltw = forcing['seaice_melt'].isel(time=n)

# liquid runoff entering ocean (non-negative)
lrunoff = forcing['lrunoff'].isel(time=n)

# frozen runoff entering ocean (non-negative)
frunoff = forcing['frunoff'].isel(time=n)

# liquid runoff from GLC entering ocean (non-negative)
lrunoff_glc = forcing['lrunoff_glc'].isel(time=n)

# frozen runoff from GLC entering ocean (non-negative)
frunoff_glc = forcing['frunoff_glc'].isel(time=n)

# liquid precipitation entering ocean.
# note: includes exchanges with sea-ice, with 
#       melt adding mass to ocean; formation removing mass.
lprec = forcing['lprec'].isel(time=n)

# frozen precipitation entering ocean.
fprec = forcing['fprec'].isel(time=n)

# virtual precipitation arising from conversion of salt restoring to water flux
vprec = forcing['vprec'].isel(time=n)

# net mass flux crossing surface (including exchange with sea-ice)
PRCmE = forcing['PRCmE'].isel(time=n)


#--------------------------------------------------------------
# heat flux crossing ocean surface and bottom [Watt/m^2]
# positive values indicate heat entering ocean;
# negative values indicate heat leaving ocean.

# geothermal heating at ocean bottom 
geothermal = forcing['internal_heat'].isel(time=n)

# net heat crossing ocean surface due to all processes, except restoring 
net_heat_surface = forcing['net_heat_surface'].isel(time=n)

# net heat passed through coupler from shortwave, longwave, latent, sensible.
# note: latent includes heat to vaporize liquid and heat to melt ice/snow. 
# note: sensible includes air-sea and ice-sea sensible heat fluxes.  
net_heat_coupler = forcing['net_heat_coupler'].isel(time=n)

# sum of longwave + latent + sensible
LwLatSens = forcing['LwLatSens'].isel(time=n)

# net shortwave passing through ocean surface
SW = forcing['SW'].isel(time=n)

# heating of liquid seawater due to formation of frazil sea ice
frazil = forcing['frazil'].isel(time=n)

# heat from melting/forming sea ice
melth = forcing['seaice_melt_heat'].isel(time=n)

# there is no restoring heat flux with this test case 
heat_restore = 0.0 * forcing['frazil'].isel(time=n)

# net heat content associated with transfer of mass across ocean surface, 
# computed relative to 0C. Both diagnostics should be the same, though 
# they are computed differently in MOM6.   
heat_pme = forcing['Heat_PmE'].isel(time=n)
heat_content_surfwater = forcing['heat_content_surfwater'].isel(time=n)

# heat content associated with water mass entering ocean
heat_content_massin = forcing['heat_content_massin'].isel(time=n)

# heat content associated with water mass leaving ocean
heat_content_massout = forcing['heat_content_massout'].isel(time=n)

# heat content associated with evaporation
heat_content_evap = forcing['heat_content_evap'].isel(time=n)

# heat content associated with meltw
# no longer needed (GMM)
#heat_content_meltw = forcing['heat_content_icemelt'].isel(time=n)

# heat content associated with liquid precipitation 
heat_content_lprec = forcing['heat_content_lprec'].isel(time=n)

# heat content associated with frozen precipitation 
heat_content_fprec = forcing['heat_content_fprec'].isel(time=n)

# heat content associated with virtual precipitation 
heat_content_vprec = forcing['heat_content_vprec'].isel(time=n)

# heat content associated with liquid runoff 
heat_content_lrunoff = forcing['heat_content_lrunoff'].isel(time=n)

# heat content associated with frozen runoff 
heat_content_frunoff = forcing['heat_content_frunoff'].isel(time=n)

# heat content associated with frozen runoff from GLC
heat_content_frunoff_glc = forcing['heat_content_frunoff_glc'].isel(time=n)

# heat content associated with liquid runoff from GLC
heat_content_lrunoff_glc = forcing['heat_content_lrunoff_glc'].isel(time=n)

# heat content associated with liquid condensation 
heat_content_cond = forcing['heat_content_cond'].isel(time=n)

#--------------------------------------------------------------
# salt flux crossing ocean surface and bottom [kg/(m^2 s)]
# positive values indicate salt entering ocean; 
# negative values indicate salt leaving ocean.

salt_flux = forcing['salt_flux_in'].isel(time=n)

# salt flux associated with surface restoring.
# salt_flux has contribution from sea ice + restoring, so we need to remove salt_flux (salt_flux_in)
salt_restore = forcing['salt_flux'].isel(time=n) - salt_flux

# for easy setup of subplots
def newSP(y, x):
    global __spv, __spi
    __spv = (y, x)
    __spi = 1
    plt.subplot(__spv[0], __spv[1], __spi)

def nextSP():
    global __spv, __spi
    __spi = __spi + 1
    plt.subplot(__spv[0], __spv[1], __spi)

def make_plot(lon, lat, field, title, xmin=-280, xmax=70, ymin=-80, ymax=90, cmin=-200, cmax=200, xlabel=False):
    '''
    Uses xarray to compute statistics and plot field as a function of lat and lon.
    Writes the max, min and ave values of field into the plot.
    '''
    global area

    # Convert to xarray DataArray if needed
    if not isinstance(field, xr.DataArray):
        field = xr.DataArray(field)

    # Compute statistics using xarray
    field_min = float(field.min())
    field_max = float(field.max())
    field_ave = float((field * area).sum() / area.sum())

    # Extract coordinate arrays
    lon_vals = lon.values if isinstance(lon, xr.DataArray) else np.asarray(lon)
    lat_vals = lat.values if isinstance(lat, xr.DataArray) else np.asarray(lat)

    # Build DataArray with 2D lon/lat as coordinates for xarray plotting
    plot_da = xr.DataArray(
        field.values,
        dims=['yh', 'xh'],
        coords={'geolon': (['yh', 'xh'], lon_vals),
                'geolat': (['yh', 'xh'], lat_vals)}
    )

    ax = plt.gca()
    kwargs = dict(x='geolon', y='geolat', ax=ax, extend='both', add_colorbar=True)
    if cmin != 0.0 or cmax != 0.0:
        kwargs['vmin'] = cmin
        kwargs['vmax'] = cmax

    plot_da.plot.pcolormesh(**kwargs)

    ax.set_title(r'' + title)
    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax)
    ax.set_ylabel(r'Latitude [$\degree$N]')
    if xlabel:
        ax.set_xlabel(r'Longitude')
    else:
        plt.xlabel('')

    ax.annotate('max=%5.2f\nmin=%5.2f\nave=%5.2f' % (field_max, field_min, field_ave),
                xy=(0.01, 0.74), xycoords='axes fraction',
                verticalalignment='bottom', fontsize=8, color='black')

---

Mass fluxes and global seawater mass budget

Global seawater mass budget consistency check

We compute the change in seawater mass over a given time period using two methods. Note that the difference from the offline diagnostics is generally larger than that from the online diagnostics, because the model output interval is usually larger than the model time step, and the loss of precision may occur depending on how the output is written by the model, especially if the output is written in single precision.

1. The net mass per time of water (units of kg/s) entering through the ocean boundaries is given by the area integral $\(\begin{equation*} \mbox{boundary water mass entering liquid seawater} = \int Q_{W} \, dA, \end{equation*}\)\( where the net water flux (units of \)\mbox{kg}~\mbox{m}^{-2}~\mbox{s}^{-1}\() is given by \)\(\begin{align*} Q_{W} &= {\tt PRCmE}. \end{align*}\)$ A nonzero surface mass flux is associated with liquid and solid precipitation and runoff; evaporation and condensation; sea ice melt/formation; and surface restoring.

2. The time change of liquid seawater mass is computed according to $\(\begin{equation*} \mbox{seawater mass change at step \)n\(} = \frac{1}{\tau_{n} - \tau_{n-1} } \int \left[\int (\rho_{n} - \rho_{n-1}) \, \mathrm{d}z \right] dA, \end{equation*}\)\( where \)\tau_{n} - \tau_{n-1}\( is the time increment of the output (i.e., one day) in seconds. Note that we make use of the MOM6 diagnostic for depth integrated density \)\(\begin{equation*} {\tt mass\_wt} = \int \rho \, \mathrm{d}z. \end{equation*}\)\( For a Boussinesq fluid, the in-situ \)\rho\( factor is set to \)\rho_{0}\(, in which case the diagnostic field \){\tt mass_wt}\( measures the thickness of a fluid column, multiplied by \)\rho_{0}\(. Note that the numerical difference of \){\tt mass_wt}\( is taken between \)\tau_{n}\( and \)\tau_{n-1}\(, because \){\tt mass_wt[n]}\( is a snapshot at the end of the \)n\(-th day, whereas \){\tt PRCmE[n]}\( is the daily average over the \)n\(-th day. For self-consistency, we should have the following equality \)\(\begin{equation*} \boxed{ \mbox{boundary water mass entering liquid seawater} = \mbox{seawater mass change}. } \end{equation*}\)$

n0       = n - 1
dmass_wt = mass_wt.isel(time=n) - mass_wt.isel(time=n0)
dt       = time[n] - time[n0]
lhs      = area * dmass_wt / dt
rhs      = area * PRCmE

# ADJUST_NET_FRESH_WATER_TO_ZERO = False

print('Total seawater mass at time step n [kg seawater]   =', float((mass_wt.isel(time=n) * area).sum()))
print('Total seawater mass at time step n0 [kg seawater]  =', float((mass_wt.isel(time=n0) * area).sum()))
print('Total seawater mass change [kg seawater]           =', float(dt * lhs.sum()))
print('Net mass flux across ocean surface [kg seawater]   =', float(dt * rhs.sum()))
print('Residual [kg seawater]                             =', float(dt * lhs.sum() - dt * rhs.sum()))
print('Residual relative to seawater mass change          =', float((lhs.sum() - rhs.sum()) / lhs.sum()))
print('Residual relative to mass_wt[n]                    =', float((lhs.sum() - rhs.sum()) / (mass_wt.isel(time=n) * area).sum()))

Total seawater mass at time step n [kg seawater]   = 1.3839157468285895e+21
Total seawater mass at time step n0 [kg seawater]  = 1.3839157948509842e+21
Total seawater mass change [kg seawater]           = -48022394659399.57
Net mass flux across ocean surface [kg seawater]   = -48022394659597.41
Residual [kg seawater]                             = 197.8359375
Residual relative to seawater mass change          = -4.119667607458938e-12
Residual relative to mass_wt[n]                    = 1.6545602858690613e-24

Surface mass fluxes I: combined fields

We will check if the net mass flux crossing surface is equal to the summation of mass flux entering and leaving ocean. If the budget is closed, then the difference should be zero.

plt.figure(figsize=(18, 8))
newSP(2, 2)

field = 86400.0 * net_massin
make_plot(lon, lat, field, 'Net mass flux entering ocean [$kg/m^2/day$]', cmin=-20, cmax=20)

nextSP()
field = 86400.0 * net_massout
make_plot(lon, lat, field, 'Net mass flux leaving ocean [$kg/m^2/day$]', cmin=-20, cmax=20)

nextSP()
field = 86400.0 * PRCmE
make_plot(lon, lat, field, 'Net mass flux across ocean surface [$kg/m^2/day$]', cmin=-20, cmax=20, xlabel=True)

nextSP()
field = 86400.0 * (PRCmE - net_massout - net_massin)
make_plot(lon, lat, field, 'Mass flux budget [$kg/m^2/day$]', cmin=-0, cmax=0, xlabel=True)

Surface mass fluxes II: component fields

plt.figure(figsize=(27, 12))
newSP(3, 3)

# lprec contains a contribution from sea ice melt/formation.
field = 86400.0 * lprec
make_plot(lon, lat, field, 'Liquid precipitation entering ocean [$kg/m^2/day$]', cmin=-20, cmax=20)

nextSP()
field = 86400.0 * fprec
make_plot(lon, lat, field, 'Frozen precipitation entering ocean [$kg/m^2/day$]', cmin=-2, cmax=2)

nextSP()
# evaporation and condensation
field = 86400.0 * evap
make_plot(lon, lat, field, 'Evaporation (negative) and condensation (positive) [$kg/m^2/day$]', cmin=-20, cmax=20)

nextSP()
field = 86400.0 * lrunoff
make_plot(lon, lat, field, 'Liquid runoff entering ocean [$kg/m^2/day$]', cmin=-2, cmax=2)

nextSP()
field = 86400.0 * frunoff
make_plot(lon, lat, field, 'Frozen runoff entering ocean [$kg/m^2/day$]', cmin=-2, cmax=2)

nextSP()
field = 86400.0 * meltw
make_plot(lon, lat, field, 'Sea ice formation (negative) and melting (positive) [$kg/m^2/day$]', cmin=-20, cmax=20)

nextSP()
field = 86400.0 * lrunoff_glc
make_plot(lon, lat, field, 'Liquid runoff from GLC entering ocean [$kg/m^2/day$]', cmin=-1, cmax=1, xlabel=True)

nextSP()
field = 86400.0 * frunoff_glc
make_plot(lon, lat, field, 'Frozen runoff from GLC entering ocean [$kg/m^2/day$]', cmin=-1, cmax=1, xlabel=True)

nextSP()
field = 86400.0 * vprec
make_plot(lon, lat, field, 'Virtual precipitation from conversion of salt restoring to water flux [$kg/m^2/day$]', cmin=-1, cmax=1, xlabel=True)

Surface mass flux self-consistency check

We will now check if the net mass flux crossing surface is equal to the summation of the above quantities. If the budget is closed, then the difference should be zero.

plt.figure(figsize=(9, 4))
newSP(1, 1)
field = 86400.0 * (PRCmE - lprec - fprec - evap - lrunoff - frunoff - meltw - frunoff_glc - lrunoff_glc - vprec)
make_plot(lon, lat, field, 'Mass flux budget [$kg/m^2/day$]', cmin=-0.0, cmax=0.0, xlabel=True)

---

Heat fluxes and global ocean heat budget

Global heat budget consistency check

We compute the change in seawater heat content over a given time period using two methods. If the difference is large, then there is a bug.

1. The net heat per time (units of Watts) entering through the ocean (surface and bottom) boundaries is given by the area integral $\(\begin{equation*} \mbox{boundary heating of liquid seawater} = \int Q \, dA, \end{equation*}\)\( where the net heat flux (units of \)\mbox{W}~\mbox{m}^{-2}~\mbox{s}^{-1}\() is given by \)\(\begin{align*} Q &= {\tt (net\_heat\_coupler + heat\_pme + frazil) + internal\_heat + heat\_restore} \\ &= {\tt net\_heat\_surface + internal\_heat + heat\_restore} \end{align*}\)\( Note that \){\tt internal_heat}$ represents the geothermal heating from the seafloor, which also contributes to the heat budget.

2. The time change of liquid seawater heat is computed according to $\(\begin{equation*} \mbox{seawater heat content change} = \frac{C_p }{\tau_{n} - \tau_{n-1} } \int dA \left(\rho_0 \int (\Theta_{n} - \Theta_{n-1}) \, \mathrm{d}z \right) \end{equation*}\)\( where \)\tau_{n} - \tau_{n-1}\( is the time increment of the output (i.e., one day) in seconds. Note that we make use of the MOM6 diagnostic for depth integrated potential/conservative temperature \)\(\begin{equation*} {\tt tomint} = \rho_0 \int \Theta \, \mathrm{d}z, \end{equation*}\)\( where the Boussinesq reference density, \)\rho_{0}\(, is used since this test case makes the Boussinesq approximation. For self-consistency, we should have the following equality \)\(\begin{equation*} \boxed{ \mbox{boundary heating of liquid seawater} = \mbox{seawater heat content change}. } \end{equation*}\)$

n0      = n - 1
dtomint = tomint.isel(time=n) - tomint.isel(time=n0)
dt      = time[n] - time[n0]
lhs     = cp * area * dtomint / dt
rhs     = area * (net_heat_coupler + heat_pme + frazil + geothermal)

print('Total seawater heat at time step n [Joules]    =', float(cp * (area * tomint.isel(time=n)).sum()))
print('Total seawater heat at time step n0 [Joules]   =', float(cp * (area * tomint.isel(time=n0)).sum()))
print('Total seawater heat content change [Joules]    =', float(dt * lhs.sum()))
print('Net heat flux across ocean boundaries [Joules] =', float(dt * rhs.sum()))
print('Residual [Joules]                              =', float(dt * lhs.sum() - dt * rhs.sum()))
print('Residual relative to heat content change       =', float((lhs.sum() - rhs.sum()) / lhs.sum()))
print('Residual relative to tomint[n]                 =', float((lhs.sum() - rhs.sum()) / (cp * (area * tomint.isel(time=n)).sum())))

Total seawater heat at time step n [Joules]    = 1.9845832072498946e+25
Total seawater heat at time step n0 [Joules]   = 1.9845176263322942e+25
Total seawater heat content change [Joules]    = 6.558091759992569e+20
Net heat flux across ocean boundaries [Joules] = 6.558091760000353e+20
Residual [Joules]                              = -778436608.0
Residual relative to heat content change       = -1.1870282217595598e-12
Residual relative to tomint[n]                 = -4.539996089398271e-22

Basic components to surface heat flux

Self-consistency check: \({\tt net\_heat\_surface = net\_heat\_coupler + heat\_pme + frazil}\). If the budget is closed, then the difference should be zero.

plt.figure(figsize=(18, 12))
newSP(3, 2)

field = net_heat_surface
make_plot(lon, lat, field, 'Net heat flux across ocean surface [$W/m^2$]')

nextSP()
field = net_heat_coupler
make_plot(lon, lat, field, 'Net heat passing through coupler [$W/m^2$]')

nextSP()
field = heat_pme
make_plot(lon, lat, field, 'Heat content associated with mass flux across ocean surface [$W/m^2$]', cmin=-20, cmax=20)

nextSP()
field = frazil
make_plot(lon, lat, field, 'Heating of liquid seawater due to frazil ice formation [$W/m^2$]', cmin=-20, cmax=20)

nextSP()
field = net_heat_surface - net_heat_coupler - heat_pme - frazil
make_plot(lon, lat, field, 'Heat flux budget [$W/m^2$]', cmin=-0.0, cmax=0.0, xlabel=True)

nextSP()
field = geothermal
make_plot(lon, lat, field, 'Geothermal heating [$W/m^2$]', cmin=-1.0, cmax=1.0, xlabel=True)

Heat fluxes acr ocean surface via the coupler

Heat flux across ocean surface via the coupler: \({\tt net\_heat\_coupler = LwLatSens + SW + melth}\), where \({\tt LwLatSens = LW + latent + sensible}\). If the budget is closed, then the difference should be zero.

plt.figure(figsize=(18, 12))
newSP(3, 2)

field = net_heat_coupler
make_plot(lon, lat, field, 'Net heat passing through coupler [$W/m^2$]', cmin=-200, cmax=200)

nextSP()
field = LwLatSens
make_plot(lon, lat, field, 'Sum of longwave, latent, and sensible heat [$W/m^2$]', cmin=-200, cmax=200)

nextSP()
field = SW
make_plot(lon, lat, field, 'Net shortwave passing through ocean surface [$W/m^2$]', cmin=-200, cmax=200)

nextSP()
field = melth
make_plot(lon, lat, field, 'Sea ice formation (positive) and melting (negative) [$W/m^2$]', cmin=-200, cmax=200, xlabel=True)

nextSP()
field = net_heat_coupler - LwLatSens - SW - melth
make_plot(lon, lat, field, 'Heat flux budget [$W/m^2$]', cmin=-0.0, cmax=0.0, xlabel=True)

Components of heat content associated with mass flux across ocean surface

Components of heat content of surface mass fluxes: $\(\begin{align*} {\tt heat\_PmE = } &{\tt heat\_content\_lprec + heat\_content\_fprec + heat\_content\_evap + heat\_content\_lrunoff + heat\_content\_frunoff +} \\ &{\tt heat\_content\_cond + heat\_content\_frunoff\_glc + heat\_content\_lrunoff\_glc + heat\_content\_vprec} \end{align*}\)$

plt.figure(figsize=(27, 12))
newSP(3, 3)

field = heat_content_lprec
make_plot(lon, lat, field, 'Heat content associated with liquid precipitation [$W/m^2$]', cmin=-20.0, cmax=20.0)

nextSP()
field = heat_content_fprec
make_plot(lon, lat, field, 'Heat content associated with frozen precipitation [$W/m^2$]', cmin=-0.2, cmax=0.2)

nextSP()
field = heat_content_evap
make_plot(lon, lat, field, 'Heat content associated with evaporation [$W/m^2$]', cmin=-20.0, cmax=20.0)

nextSP()
field = heat_content_lrunoff
make_plot(lon, lat, field, 'Heat content associated with liquid runoff [$W/m^2$]', cmin=-1.0, cmax=1.0)

nextSP()
field = heat_content_frunoff
make_plot(lon, lat, field, 'Heat content associated with frozen runoff [$W/m^2$]', cmin=-0.2, cmax=0.2)

nextSP()
field = heat_content_cond
make_plot(lon, lat, field, 'Heat content associated with liquid condensation [$W/m^2$]', cmin=-1.0, cmax=1.0)

nextSP()
field = heat_content_lrunoff_glc
make_plot(lon, lat, field, 'Heat content associated with liquid runoff from GLC [$W/m^2$]', cmin=-0.0, cmax=0.0, xlabel=True)

nextSP()
field = heat_content_frunoff_glc
make_plot(lon, lat, field, 'Heat content associated with frozen runoff from GLC [$W/m^2$]', cmin=-0.0, cmax=0.0, xlabel=True)

nextSP()
field = heat_content_vprec
make_plot(lon, lat, field, 'Heat content associated with virtual precipitation [$W/m^2$]', cmin=-0.0, cmax=0.0, xlabel=True)

Self-consistency of heat content associated with mass flux across ocean surface

comp_sum = (heat_content_lprec + heat_content_fprec + heat_content_evap +
            heat_content_lrunoff + heat_content_frunoff + heat_content_cond +
            heat_content_lrunoff_glc + heat_content_frunoff_glc + heat_content_vprec)

plt.figure(figsize=(18, 12))
newSP(3, 2)

field = heat_content_surfwater
make_plot(lon, lat, field, '(a) heat content associated with mass flux across ocean surface [$W/m^2$]', cmin=-20.0, cmax=20.0)

nextSP()
field = heat_pme
make_plot(lon, lat, field, '(b) Heat flux into ocean from mass flux into ocean [$W/m^2$]', cmin=-20.0, cmax=20.0)

nextSP()
field = heat_pme - heat_content_surfwater
make_plot(lon, lat, field, '(c) Difference between (b) and (a) [$W/m^2$]', cmin=0.0, cmax=0.0)

nextSP()
field = comp_sum
make_plot(lon, lat, field, '(d) Sum of heat content components [$W/m^2$]', cmin=-20.0, cmax=20.0, xlabel=True)

nextSP()
field = comp_sum - heat_content_surfwater
make_plot(lon, lat, field, '(e) Difference between (d) and (a) [$W/m^2$]', cmin=0.0, cmax=0.0, xlabel=True)

---

Salt fluxes and global ocean salt budget

Global salt budget consistency check

We compute the change in seawater salt content over a given time period using two methods. If the difference is large, then there is a bug.

1. The net salt per unit time (units of kg/s) entering through the ocean boundaries is given by the area integral $\(\begin{equation*} \mbox{boundary salt entering liquid seawater} = \int Q_{S} \, dA, \end{equation*}\)\( where the net salt flux (units of \)\mbox{kg}~\mbox{m}^{-2}~\mbox{s}^{-1}\() is given by \)\(\begin{align*} Q_{S} &= {\tt salt\_flux} + {\tt salt\_restore}. \end{align*}\)$ A nonzero salt flux is associated with exchanges between liquid seawater and solid sea ice. It also arises from simulations using a restoring boundary flux associated with damping to observed sea surface salinity. Finally, there can be a salt flux when using sponges to damp the ocean interior back to an observed value.

2. The time change of liquid seawater salt content is computed according to $\(\begin{equation*} \mbox{seawater salt content change} = \frac{1}{\tau_{n} - \tau_{n-1} } \int dA \left(\rho_0 \int (S_{n} - S_{n-1}) \, \mathrm{d}z \right) \end{equation*}\)\( where \)\tau_{n} - \tau_{n-1}\( is the time increment of the output (i.e., one day) in seconds. Note that we make use of the MOM6 diagnostic for depth integrated salinity \)\(\begin{equation*} {\tt somint} = \rho_0 \int S \, \mathrm{d}z, \end{equation*}\)\( where the Boussinesq reference density, \)\rho_{0}\(, is used since this test case makes the Boussinesq approximation. For self-consistency, we should have the following equality \)\(\begin{equation*} \boxed{ \mbox{boundary salt entering liquid seawater} = \mbox{seawater salt content change}. } \end{equation*}\)$

n0      = n - 1
dsomint = somint.isel(time=n) - somint.isel(time=n0)
dt      = time[n] - time[n0]
lhs     = 1.e-3 * area * dsomint / dt
rhs     = area * (salt_flux + salt_restore)

print('Total seawater salt at time step n [kg salt]     =', float((area * somint.isel(time=n)).sum()))
print('Total seawater salt at time step n0 [kg salt]    =', float((area * somint.isel(time=n0)).sum()))
print('Total seawater salt content change [kg salt]     =', float(dt * lhs.sum()))
print('Net salt flux across ocean surface [kg salt]     =', float(dt * rhs.sum()))
print('Residual [kg salt]                               =', float(dt * lhs.sum() - dt * rhs.sum()))
print('Residual relative to salt content change         =', float((lhs.sum() - rhs.sum()) / lhs.sum()))
print('Residual relative to somint[n]                   =', float((lhs.sum() - rhs.sum()) / (area * somint.isel(time=n)).sum()))

Total seawater salt at time step n [kg salt]     = 4.805000740653248e+22
Total seawater salt at time step n0 [kg salt]    = 4.805000745001683e+22
Total seawater salt content change [kg salt]     = -43484349186.481476
Net salt flux across ocean surface [kg salt]     = -43484349598.13806
Residual [kg salt]                               = 411.6565856933594
Residual relative to salt content change         = -9.466775886436414e-09
Residual relative to somint[n]                   = 9.915802522023821e-26

plt.figure(figsize=(18, 8))
newSP(2, 2)

field = 86400.0 * salt_flux
make_plot(lon, lat, field, 'Salt flux crossing ocean surface and bottom [kg m$^{-2}$ day$^{-1}$]', cmin=-1.0, cmax=1.0, xlabel=True)

nextSP()
field = 86400.0 * salt_restore
make_plot(lon, lat, field, 'Salt flux associated with surface restoring [kg m$^{-2}$ day$^{-1}$]', cmin=-0.1, cmax=0.1, xlabel=True)