Compute column inventories of tritium at BATS

Ivan Lima - Thu Mar 29 2018 10:56:28 -0400

In [1]:
%matplotlib notebook
import pandas as pd
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import os, gsw
from matplotlib import colors
from cesm_utils import da2ma, fixtime, spd, dpm
from bgc_utils import rho
from scipy import interpolate
pd.options.display.max_columns = 50

Read database

In [2]:
def dyear2date(strn_list, yr0=1950): # convert decimal year to date objects
    ndays = [np.round((np.float(x)-yr0) * 365.2425) for x in strn_list] # number of days since 1950-01-01
    return pd.to_datetime(ndays, unit='D',origin='%d-01-01'%yr0)

database = pd.read_csv('data/tritium-helium_database/tritium_helium_all.csv', index_col=0, na_values=-999,
                       parse_dates=[6], date_parser=dyear2date, low_memory=False)
database['depth'] = gsw.z_from_p(np.abs(database.pressure), database.latitude) # compute depth from pressure

Funtion declarations

In [3]:
def get_box(df, lon0, lat0, delta=1):
    """ get data in box centered at lon0, lat0"""
    box = df[df.latitude < lat0+delta]
    box = box[box.latitude > lat0-delta]
    box = box[box.longitude < lon0+delta]
    box = box[box.longitude > lon0-delta]
    return box.sort_values('stadate') # sort by date

# depth bins
zbins = np.array([-4500, -4000, -3500, -3000, -2500, -2000, -1750, -1500, -1250, -1000, -900, -800, -700,
                  -600, -500, -400, -300, -250, -200, -150, -100, -50, 0])

def get_profile(df, vname='tritium', xvar='stadate'):
    """ bin-average by depth to create series of vertical profiles (depth x time) """
    factor = pd.cut(df.depth, bins=zbins)
    df_binned = df.groupby([factor,xvar])[vname].mean().unstack()
    if vname in ['tritium']:
        df_binned = df_binned * 0.10995 # TU -> pmol/m^3
    df_binned = df_binned.resample('M', axis=1).mean()
    x = df_binned.columns
    return df_binned, x.insert(524,x[-1] + pd.Timedelta(30,'D')), zbins

def get_inventory(df, nn=17):
    """ compute time series of tritium inventory (pmol/m^2)"""
    df_work = df.copy()
    df_work.loc[:,df_work.count() < nn] = np.nan # stations with less that 17 vertical points are ignored
#     dz = np.diff(zbins[5:])                    # integrate down to 2000 m
#     inventory = df_work.iloc[5:,:].mul(dz,axis='index').sum(axis=0)
    dz = np.diff(zbins)                    # integrate down to 2000 m
    inventory = df_work.mul(dz,axis='index').sum(axis=0)
    inventory.name = 'observations'
    inventory[inventory==0] = np.nan
    return inventory.sort_index() # sort by time

Ir = 1.384e-6 # 3He/4He isotopic ratio

from tr3he import alpha_sol_he, comp_he_sol_0, comp_sol_0

def get_excess_helium3(temp,sal,helium3,neon):
    alpha_sol = alpha_sol_he(temp)
    he4_ceq = comp_he_sol_0(temp,sal,rho) / rho * 1.e+9   # 4He equilibrium concentration (mol/m^3 -> nmol/kg)
    ne_ceq  = comp_sol_0(temp,sal,rho,'Ne') / rho * 1.e+9 # Ne equilibrium concentration (mol/m^3 -> nmol/kg)
    he3_ceq = he4_ceq * Ir  * alpha_sol                   # 3He equilibrium concentration
    if neon is not None:
        Ca_he4 = 0.28823 * (neon - ne_ceq)
    else:
        Ca_he4 = 0.003156 * ne_ceq
    Ca_he3 = Ir * Ca_he4
    return helium3 - he3_ceq - Ca_he3

Extract observational data at BATS location & compute tritium inventory

In [4]:
bats_lat = database[database.sect_id=='BATS'].latitude.mean()
bats_lon = database[database.sect_id=='BATS'].longitude.mean()
bats = get_box(database, bats_lon, bats_lat, delta=2) # all data in a 2x2 degree box
bats_profile, x, y = get_profile(bats)
bats_inventory = get_inventory(bats_profile)

Read model output at BATS & compute tritium inventory

In [5]:
ds_bats = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/BATS.tr3he.GIAF.01.nc', decode_times=False)
ds_bats = fixtime(ds_bats, 1824)
ds_bats = ds_bats.squeeze(dim=('nlat','nlon'))
ds_bats = ds_bats.drop(['TLONG','TLAT'])
model_tritium_inv = (ds_bats.TRITIUM * ds_bats.dz * 1.e-2).sum(dim='z_t').to_series() # cm -> m

Plot figures

Tritium profiles & inventory from observations

Profiles with less than 17 depth bins are discarded in the calculation of the inventory.

In [6]:
fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(hspace=0.05,right=0.875)
_ = axs[0].plot(bats_inventory.index, bats_inventory.values, '.', label='observations')
_ = axs[0].set(ylabel='tritium inventory (pmol m$^{-2}$)', title='BATS')
xmin, xmax = axs[0].get_xlim()

pm = axs[1].pcolormesh(x, y, bats_profile, cmap=plt.cm.plasma, norm=colors.LogNorm())
_ = axs[1].set(xlabel='time', ylabel='depth (m)', xlim=(xmin, xmax))
l, b, w, h =  axs[1].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritium concentration (pmol m$^{-3}$)')
# _ = axs[1].plot(bats.stadate, bats.depth, 'x')

Tritium inventory from model and observations

In [7]:
fig, ax = plt.subplots()
_ = bats_inventory.plot(ax=ax, style='.', label='observations')
# _ = bats_inventory.resample('A').mean().plot(ax=ax, style='.', label='observations')
_ = model_tritium_inv.plot(ax=ax, style='.', label='model')
_ = ax.set(ylabel='tritium inventory (pmol m$^{-2}$)', title='BATS')
_ = ax.legend(loc='best')

Surface tritium from model and observations

In [8]:
bats_surface = bats_profile.iloc[-1,:]
model_tritium_surf = ds_bats.TRITIUM.sel(z_t=500).to_series()
fig, ax = plt.subplots()
_ = bats_surface.plot(ax=ax, style='.', label='observations')
_ = model_tritium_surf.plot(ax=ax, style='.', label='model')
_ = ax.set(ylabel='surface tritium (pmol m$^{-3}$)', title='BATS')
_ = ax.legend(loc='best')

Tritium profiles from model and observations

In [9]:
fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(right=0.875)
pm = axs[0].pcolormesh(x, y, bats_profile, cmap=plt.cm.plasma, norm=colors.LogNorm())
_ = axs[0].set(ylabel='depth (m)', title='observations')
l, b, w, h =  axs[0].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritium (pmol m$^{-3}$)')
ymin, ymax = axs[0].get_ylim()
vmin, vmax = pm.get_clim()
data = da2ma(ds_bats.TRITIUM.transpose().values)
data[data<=0] = vmin
pm = axs[1].pcolormesh(ds_bats.time, -ds_bats.z_t/100, data, cmap=plt.cm.plasma, vmin=vmin, vmax=vmax, 
                       norm=colors.LogNorm())
_ = axs[1].set(xlabel='time', ylabel='depth (m)', ylim=(ymin,ymax), title='model')

Temperature, salinity, helium and helium 3 profiles from observations

Those are the input variables for the excess helium 3 calculation.

In [10]:
delhe3, x, y = get_profile(bats,'delhe3')
helium, x, y = get_profile(bats,'helium')
temp,   x, y = get_profile(bats,'temperature')
salt,   x, y = get_profile(bats,'salinity')
helium3 = Ir * (1 + delhe3/100) * helium

import cmocean
fig, axs = plt.subplots(2,2,sharex=True, sharey=True, figsize=(9.5,7.2))
fig.subplots_adjust(wspace=0.075)
pm = axs[0,0].pcolormesh(x, y, temp, cmap=cmocean.cm.thermal)
cb = fig.colorbar(pm,ax=axs[0,0])
cb.set_label(r'$^{\circ}$C')
_ = axs[0,0].set_title('temperature')

pm = axs[0,1].pcolormesh(x, y, salt, cmap=cmocean.cm.haline)
cb = fig.colorbar(pm,ax=axs[0,1])
cb.set_label(r'g/kg')
_ = axs[0,1].set_title('salinity')

pm = axs[1,0].pcolormesh(x, y, helium, cmap=plt.cm.plasma)
cb = fig.colorbar(pm,ax=axs[1,0])
cb.set_label('nmol kg$^{-1}$')
_ = axs[1,0].set_title('helium')

pm = axs[1,1].pcolormesh(x, y, helium3, cmap=plt.cm.plasma)
cb = fig.colorbar(pm,ax=axs[1,1])
cb.set_label('nmol kg$^{-1}$')
_ = axs[1,1].set_title('helium 3')

for i in range(2):
    _ = axs[i,0].set_ylabel('depth (m)')

for j in range(2):
    plt.setp(axs[1,j].get_xticklabels(),rotation=35,horizontalalignment='right')

Tritiumgenic helium 3 profiles & inventory from observations

Profiles with less than 17 depth bins are discarded in the calculation of the inventory.

In [11]:
excess_helium3 = get_excess_helium3(temp, salt, helium3, None) * rho * 1.e+3 # nmol/kg -> pmol/m^3
excess_helium3[excess_helium3<0] = np.nan
excess_helium3_inv = get_inventory(excess_helium3)

fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(hspace=0.05,right=0.875)
_ = axs[0].plot(excess_helium3_inv.index, excess_helium3_inv.values, '.', label='observations')
_ = axs[0].set(ylabel='tritiumgenic helium 3 inventory (pmol m$^{-2}$)', title='BATS')
xmin, xmax = axs[0].get_xlim()

pm = axs[1].pcolormesh(x, y, excess_helium3, cmap=plt.cm.plasma)
_ = axs[1].set(xlabel='time', ylabel='depth (m)', xlim=(xmin, xmax))
l, b, w, h =  axs[1].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritiumgenic helium 3 (pmol m$^{-3}$)')

Tritiumgenic helium 3 profiles from model and observations

In [12]:
model_excess_helium3 = ds_bats.HELIUM3 - ds_bats.INERT_HELIUM3

fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(right=0.875)
pm = axs[0].pcolormesh(x, y, excess_helium3, cmap=plt.cm.plasma, norm=colors.LogNorm())
_ = axs[0].set(ylabel='depth (m)', title='observations')
l, b, w, h =  axs[0].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritiumgenic helium 3 (pmol m$^{-3}$)')
ymin, ymax = axs[0].get_ylim()
vmin, vmax = pm.get_clim()
data = da2ma(model_excess_helium3.transpose().values)
data[data<=0] = vmin
pm = axs[1].pcolormesh(ds_bats.time, -ds_bats.z_t/100, data, cmap=plt.cm.plasma, vmin=vmin, vmax=vmax,
                       norm=colors.LogNorm())
_ = axs[1].set(xlabel='time', ylabel='depth (m)', ylim=(ymin,ymax), title='model')

Tritiumgenic helium 3 inventory from model and observations

In [13]:
model_excess_helium3_inv = (model_excess_helium3 * ds_bats.dz * 1.e-2).sum(dim='z_t').to_series() # cm -> m
fig, ax = plt.subplots()
_ = excess_helium3_inv.plot(ax=ax, style='.', label='observations')
_ = model_excess_helium3_inv.plot(ax=ax, style='.', label='model')
_ = ax.set(ylabel='tritiumgenic helium3 (pmol m$^{-2}$)', title='BATS')
_ = ax.legend(loc='best')

Model tritium budget at BATS

In [14]:
decay_vint = (ds_bats.TRITIUM_DECAY * ds_bats.dz/100 * spd).sum(dim='z_t').to_series() # pmol/m^2/s -> pmol/m^2/day
stf_tritium = (ds_bats.STF_TRITIUM * spd).to_series()                                  # pmol/m^2/s -> pmol/m^2/day
dtritium_dt = model_tritium_inv.diff()/np.resize(dpm, len(stf_tritium))
model_bats_tritium = pd.DataFrame({'stf':stf_tritium, 'decay':decay_vint, 'd3H/dt':dtritium_dt})
model_bats_tritium['stf - decay'] = model_bats_tritium.stf - model_bats_tritium.decay
model_bats_tritium['d3H/dt - (stf-decay)'] = model_bats_tritium['d3H/dt'] - model_bats_tritium['stf - decay']
model_bats_tritium = model_bats_tritium.groupby(model_bats_tritium.index.year).mean()

The monthly $\Delta^3\!H/\Delta t$ is very noisy so the plot shows annual means. STF = surface tracer flux.

In [15]:
plt.rcParams['xtick.minor.bottom'] = False # fix strange minor ticks in time series plots
fig, axs = plt.subplots(2,2,sharex=True,figsize=(9.5,7.2))
fig.subplots_adjust(hspace=0.1)
_ = model_bats_tritium.stf.plot(ax=axs[0,0], legend=True)
_ = model_bats_tritium.decay.plot(ax=axs[0,1], legend=True)
_ = model_bats_tritium[['stf - decay','d3H/dt']].plot(ax=axs[1,0], legend=True, alpha=0.75)
x, y1, y2 = model_bats_tritium.index, model_bats_tritium['d3H/dt'], model_bats_tritium['stf - decay']
_ = axs[1,0].fill_between(x, y1, y2, where=y1>=y2, facecolor='C1', alpha=0.25)
_ = axs[1,0].fill_between(x, y1, y2, where=y1<y2, facecolor='C0', alpha=0.25)
_ = model_bats_tritium['d3H/dt - (stf-decay)'].plot(ax=axs[1,1], legend=True)
_ = axs[1,0].hlines(0,model_bats_tritium.index[0], model_bats_tritium.index[-1], linestyles='dotted')
_ = axs[1,1].hlines(0,model_bats_tritium.index[0], model_bats_tritium.index[-1], linestyles='dotted')
for i in range(2):
    _ = axs[i,0].set_ylabel('pmol m$^{-2}$ d$^{-1}$')

Apply smoothing spline to observations

In [16]:
zt = (zbins[:-1] + zbins[1:])/2
ztnew = np.arange(zt.min(),0,25)
dz = np.diff(ztnew)[0]
zplot = np.arange(ztnew[0]-dz/2,0,dz)
dznew = np.diff(zplot)

def smooth_spline(df):
    profiles = np.ma.masked_where(np.isnan(df.values), df.values)
    profiles_smooth = np.ma.masked_all((len(ztnew), profiles.shape[-1]))
    for j in range(profiles.shape[-1]):
        y = np.ma.masked_where(np.isnan(profiles[:,j]), profiles[:,j]).compressed()
        x = np.ma.masked_where(np.isnan(profiles[:,j]), zt).compressed()
        s = 0.10075 - 0.00475 * len(x) 
        if len(y) > 15:
            tck = interpolate.splrep(x,y,s=s)
            profiles_smooth[:,j] = interpolate.splev(ztnew, tck, ext=3)

    return pd.DataFrame(profiles_smooth, index=ztnew, columns=df.columns)

bats_profile_smooth = smooth_spline(bats_profile)
bats_inventory_smooth = bats_profile_smooth.mul(dznew,axis='index').sum(axis=0)
bats_inventory_smooth[bats_inventory_smooth==0] = np.nan

Plot smoothed profiles and inventories

In [17]:
x = bats_profile_smooth.columns
x = x.insert(524,x[-1] + pd.Timedelta(30,'D'))
y = zplot

fig, axs = plt.subplots(2,1,sharex=True,figsize=(9.5,9))
fig.subplots_adjust(hspace=0.05,right=0.875)
_ = axs[0].plot(bats_inventory_smooth.index, bats_inventory_smooth.values, '.', label='observations')
_ = axs[0].set(ylabel='tritium inventory (pmol m$^{-2}$)', title='BATS')
xmin, xmax = axs[0].get_xlim()

pm = axs[1].pcolormesh(x, y, bats_profile_smooth, cmap=plt.cm.plasma, norm=colors.LogNorm())
_ = axs[1].set(xlabel='time', ylabel='depth (m)', xlim=(xmin, xmax))
l, b, w, h =  axs[1].get_position().bounds
cax = fig.add_axes((l+w+0.02, b, 0.02, h))
cb = fig.colorbar(pm, cax=cax)
cb.set_label('tritium concentration (pmol m$^{-3}$)')
In [18]:
# zt = (zbins[:-1] + zbins[1:])/2
# ztnew = np.arange(zt.min(),0,25)
# fig, axs = plt.subplots(1, 4, sharey=True, figsize=(9.75,4.8))
# fig.subplots_adjust(wspace=0.05)
# for date, ax in zip(['1972-06-30','1977-07-31','1988-03-31','1997-09-30'], axs.ravel()):
#     prof = bats_profile.loc[:,date].values
#     y = np.ma.masked_where(np.isnan(prof), prof).compressed()
#     x = np.ma.masked_where(np.isnan(prof), zt).compressed()
#     s = 0.10075 - 0.00475 * len(y)
#     tck = interpolate.splrep(x,y,s=s)
#     prof_new = interpolate.splev(ztnew, tck, ext=3)
#     _ = ax.plot(prof, zt, '.')
#     _ = ax.plot(prof_new, ztnew, alpha=0.5)
#     _ = ax.set(title=date)

Compute inventories for 0-1000m, 1000-3000m & >3000m

In [19]:
bats_inv_0_1000 = (bats_profile_smooth[bats_profile_smooth.index>=-1000] * dz).sum(axis=0)
bats_inv_1000_3000 = (bats_profile_smooth[(bats_profile_smooth.index<-1000) & (bats_profile_smooth.index>=-3000)]
                      * dz).sum(axis=0)
bats_inv_3000 = (bats_profile_smooth[bats_profile_smooth.index<-3000] * dz).sum(axis=0)
bats_inv_0_1000[bats_inv_0_1000==0] = np.nan
bats_inv_1000_3000[bats_inv_1000_3000==0] = np.nan
bats_inv_3000[bats_inv_3000==0] = np.nan

ds_bats['TRITIUM'] = xr.where(ds_bats.TRITIUM<0,0,ds_bats.TRITIUM)
model_inv_0_1000 = (ds_bats.TRITIUM[:,ds_bats.z_t<100000] * ds_bats.dz[ds_bats.z_t<100000]/100).sum(dim='z_t').to_series()
model_inv_1000_3000 = (ds_bats.TRITIUM[:,(ds_bats.z_t>100000) & (ds_bats.z_t<=300000)] * 
                       ds_bats.dz[(ds_bats.z_t>100000) & (ds_bats.z_t<=300000)]/100).sum(dim='z_t').to_series()
model_inv_3000 = (ds_bats.TRITIUM[:,ds_bats.z_t>300000] * ds_bats.dz[ds_bats.z_t>300000]/100).sum(dim='z_t').to_series()

fig, axs = plt.subplots(2, 2, figsize=(9.75,8))
fig.subplots_adjust(hspace=0.25)
fig.delaxes(axs[1,1])
_ = bats_inv_0_1000.plot(ax=axs[0,0], style='.', label='observations')
_ = model_inv_0_1000.plot(ax=axs[0,0], style='.', label='model')
_ = axs[0,0].set(ylabel='tritium inventory (pmol m$^{-2}$)',title='BATS 0-1000 m')
_ = axs[0,0].legend(loc='upper right')
_ = bats_inv_1000_3000.plot(ax=axs[0,1], style='.', label='observations')
_ = model_inv_1000_3000.plot(ax=axs[0,1], style='.', label='model')
_ = axs[0,1].set(title='BATS 1000-3000 m')
_ = bats_inv_3000.plot(ax=axs[1,0], style='.', label='observations')
_ = model_inv_3000.plot(ax=axs[1,0], style='.', label='model')
_ = axs[1,0].set(ylabel='tritium inventory (pmol m$^{-2}$)',title='BATS >3000 m')