Model-data comparison

Created by Ivan Lima on Mon, 05 Feb 2018 14:03:17 -0500

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 cesm_utils import da2ma
from bgc_utils import rho
import warnings
from datetime import datetime
from cases import *
print('Last updated on %s'%datetime.now().ctime())
warnings.filterwarnings('ignore')
pd.options.display.max_columns = 50

Last updated on Tue Apr 30 09:31:18 2019

Case description

In [2]:
case = 'tr3he.GIAF.06'
print(case,':',case_info[case])

tr3he.GIAF.06 : 186-year run with interannualy varying forcing (CORE 1948-2009). Physics is
spun-up for 126 years (2 IAF cycles: 1824-1949). Model is spun-up with pre-bomb tritium fluxes and
historic tritium deposition starts in the third IAF cycle (1950-2009). Time & latitude dependent
tritium concentration in precipitation (Cp) is estimated from GNIP data (Michaela's blended fits).
Tritium concentration in water vapor (Cv) is set to 0.7 Cp. Relative humidity has a max value of
0.95. Added tritium river fluxes. Time & latitude dependent tritium concentration in rivers is
estimated from GNIP data (Michaela's blended fits). Tritium concentration in precipitation over
ocean and land are smoothed over a 1 degree latitudinal grid. Run includes vapor flux correction.

Read tritium-helium3 database

In [3]:
def dyear2date(strn_list, yr0=1950):
    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

Depth ranges for binning observations & some function definitions

In [4]:
zbins = np.array([-4500, -4000, -3500, -3000, -2500, -2000, -1750, -1500, -1250, -1000, -900, -800, -700,
                  -600, -500, -450, -400, -350, -300, -250, -200, -150, -100, -75, -50, -25, 0])

info_cols = ['sect_id', 'station', 'castno', 'stadate','latitude', 'longitude','depth']
data_cols = ['tritium','tritium_flag', 'helium', 'helium_flag', 'neon', 'neon_flag']

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

def 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
In [5]:
# fig, axs = plt.subplots(1,2,figsize=(9.5,4.8))
# _ = axs[0].hist(data.depth[data.depth.notnull()],zbins=len(zbins),range=(-4500,0),rwidth=1)
# _ = axs[1].hist(data.depth[data.depth.notnull()],zbins=zbins,rwidth=1)

# np.sort(database.sect_id.unique()[1:])

# varname = 'tritium'
# data = database[(database.sect_id=='P16N') & (database.stadate.dt.year==2006) & 
#                 ((database[varname+'_flag']==2) | (database[varname+'_flag']==6))][info_cols + [varname]]
# print(data.stadate.dt.year.unique())
# print(data.stadate.min(), data.stadate.max())
# print(data.latitude.min(), data.latitude.max())
# print(data.depth.min(), data.depth.max())

BATS

Observations are bin-averaged by depth and month.

In [6]:
def get_bats_section(sec,varname,xvar='stadate'):
    if varname in ['temperature','salinity']:
        data = database[database.sect_id==sec][info_cols + [varname]]
    else:
        data = database[(database.sect_id==sec) & 
                        ((database[varname+'_flag']==2) | (database[varname+'_flag']==6))][info_cols + [varname]]
        # data = database[(database.sect_id==sec)][info_cols + [varname]]
        
    factor = pd.cut(data.depth, bins=zbins)
    # y = data.groupby(factor)['depth'].mean().values
    y = zbins
    data = data.groupby([factor,xvar])[varname].mean().unstack()
    data = data.resample('M',axis=1).mean() # bin-average by month
    x = data.columns.values
    return data, x, y

Read model results

In [7]:
model_bats = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/%s_BATS.nc'%case)
model_bats['time'] = model_bats.time.to_index().to_datetimeindex()

Tritium

In [8]:
bats_tritium, bats_dates, bats_z = get_bats_section('BATS','tritium')
model_tritium = model_bats.TRITIUM.transpose().values/0.10995 # pmol/m^3 -> TU
model_tritium[model_tritium<0] = 0

fig, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(9.5,4.8))
fig.subplots_adjust(left=0.1, right=0.9, wspace=0.05)
ymin, ymax = -4150, 0
# vmin = min(da2ma(bats_tritium.values).min(), da2ma(model_tritium).min())
# vmax = max(da2ma(bats_tritium.values).max(), da2ma(model_tritium).max())

pm = axs[0].pcolormesh(bats_dates, bats_z, bats_tritium.values[:,1:], cmap=plt.cm.plasma)#, vmin=vmin, vmax=vmax)
_ = axs[0].set(xlabel='time',ylabel='depth (m)',title='BATS tritium - observations')
plt.setp(axs[0].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[0].set_xlim(bats_dates.min(),bats_dates.max())
_ = axs[0].set_ylim(ymin,ymax)
vmin, vmax = pm.get_clim()

pm = axs[1].pcolormesh(model_bats.time, -model_bats.z_t/100, model_tritium, cmap=plt.cm.plasma, vmin=vmin, vmax=vmax)
plt.setp(axs[1].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[1].set(xlabel='time',title='BATS tritium - model')

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('TU')

Helium

In [9]:
bats_helium, bats_dates, bats_z = get_bats_section('BATS','helium')
model_helium = model_bats.HELIUM4.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_helium[model_helium<0] = 0

fig, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(9.5,4.8))
fig.subplots_adjust(left=0.1, right=0.875, wspace=0.05)
ymin, ymax = -4150, 0
# vmin = min(da2ma(bats_helium.values).min(), da2ma(model_helium).min())
# vmax = max(da2ma(bats_helium.values).max(), da2ma(model_helium).max())

pm = axs[0].pcolormesh(bats_dates, bats_z, bats_helium.values[:,1:], cmap=plt.cm.plasma)#, vmin=vmin, vmax=vmax)
_ = axs[0].set(xlabel='time',ylabel='depth (m)',title='BATS He')
plt.setp(axs[0].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[0].set_xlim(bats_dates.min(),bats_dates.max())
_ = axs[0].set_ylim(ymin,ymax)
vmin, vmax = pm.get_clim()

pm = axs[1].pcolormesh(model_bats.time, -model_bats.z_t/100, model_helium, cmap=plt.cm.plasma, vmin=vmin, vmax=vmax)
plt.setp(axs[1].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[1].set(xlabel='time',title='BATS helium - model')

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('nmol kg$^{-1}$')

Helium 3

In [10]:
bats_delhe3, bats_dates, bats_z = get_bats_section('BATS','delhe3')
bats_helium3 = Ir * (1 + bats_delhe3/100) * bats_helium
model_helium3 = model_bats.HELIUM3.transpose().values / rho * 1.e-3 # pmol/m^3 -> nmol/kg
model_helium3[model_helium3<0] = 0

fig, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(9.5,4.8))
fig.subplots_adjust(left=0.1, right=0.85, wspace=0.05)
ymin, ymax = -4150, 0
# vmin = min(da2ma(bats_helium.values).min(), da2ma(model_helium).min())
# vmax = max(da2ma(bats_helium.values).max(), da2ma(model_helium).max())

pm = axs[0].pcolormesh(bats_dates, bats_z, bats_helium3.values[:,1:], cmap=plt.cm.plasma)
_ = axs[0].set(xlabel='time',ylabel='depth (m)',title='BATS $^3$He')
plt.setp(axs[0].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[0].set_xlim(bats_dates.min(),bats_dates.max())
_ = axs[0].set_ylim(ymin,ymax)
vmin, vmax = pm.get_clim()

pm = axs[1].pcolormesh(model_bats.time, -model_bats.z_t/100, model_helium3, cmap=plt.cm.plasma, vmin=vmin, vmax=vmax)
plt.setp(axs[1].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[1].set(xlabel='time',title='BATS $^3$He - model')

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('nmol kg$^{-1}$')

Excess helium 3

In [11]:
from tr3he import alpha_sol_he, comp_he_sol_0, comp_sol_0
bats_temp, _, _ = get_bats_section('BATS','temperature')
bats_sal,  _, _ = get_bats_section('BATS','salinity')

bats_excess_helium3 = excess_helium3(bats_temp,bats_sal,bats_helium3,None)
bats_excess_helium3[bats_excess_helium3<0] = 0

model_inert_helium3 = model_bats.INERT_HELIUM3.transpose().values 
model_helium3 = model_bats.HELIUM3.transpose().values
model_excess_helium3 = (model_helium3 - model_inert_helium3) / rho * 1.e-3 # pmol/m^3 -> nmol/kg
model_excess_helium3[model_excess_helium3<0] = 0

fig, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(9.5,4.8))
fig.subplots_adjust(left=0.1, right=0.85, wspace=0.05)
ymin, ymax = -4150, 0

pm = axs[0].pcolormesh(bats_dates, bats_z, bats_excess_helium3.values[:,1:], cmap=plt.cm.plasma)
_ = axs[0].set(xlabel='time',ylabel='depth (m)',title='BATS excess $^3$He')
plt.setp(axs[0].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[0].set_xlim(bats_dates.min(),bats_dates.max())
_ = axs[0].set_ylim(ymin,ymax)
vmin, vmax = pm.get_clim()

pm = axs[1].pcolormesh(model_bats.time, -model_bats.z_t/100, model_excess_helium3, cmap=plt.cm.plasma,
                       vmin=vmin, vmax=vmax)
plt.setp(axs[1].get_xticklabels(),rotation=35,horizontalalignment='right')
_ = axs[1].set(xlabel='time',title='BATS excess $^3$He - model')

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('nmol kg$^{-1}$')

WOCE sections

A16N

Observations are bin-averaged by depth.

In [12]:
def get_woce_section(sec,year,varname,xvar):
    data = database[(database.sect_id==sec) & (database.stadate.dt.year==year) & 
                             ((database[varname+'_flag']==2) | (database[varname+'_flag']==6))][info_cols + [varname]]
    factor = pd.cut(data.depth, bins=zbins)
    # y = data.groupby(factor)['depth'].mean().values
    y = zbins
    data = data.groupby([factor,xvar])[varname].mean().unstack()
    x = data.columns.values
    return data, x, y

def get_woce_section_ts(sec,year,varname,xvar): # for TEMP & SAL only
    data = database[(database.sect_id==sec) & (database.stadate.dt.year==year) & 
                             ((database['helium_flag']==2) | (database['helium_flag']==6))][info_cols + [varname]]
    factor = pd.cut(data.depth, bins=zbins)
    # y = data.groupby(factor)['depth'].mean().values
    y = zbins
    data = data.groupby([factor,xvar])[varname].mean().unstack()
    x = data.columns.values
    return data, x, y

Read model results

In [13]:
model_A16N_1988 = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/%s_A16N_1988.nc'%case)
model_A16N_2003 = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/%s_A16N_2003.nc'%case)

Tritium

In [14]:
A16N_1988_tritium, A16N_1988_lat, A16N_1988_z = get_woce_section('A16N',1988,'tritium','latitude')
A16N_2003_tritium, A16N_2003_lat, A16N_2003_z = get_woce_section('A16N',2003,'tritium','latitude')

model_tritium_1988 = model_A16N_1988.TRITIUM.transpose().values/0.10995 # pmol/m^3 -> TU
model_tritium_2003 = model_A16N_2003.TRITIUM.transpose().values/0.10995 # pmol/m^3 -> TU
model_tritium_1988[model_tritium_1988<0] = 0
model_tritium_2003[model_tritium_2003<0] = 0

vmin = min(A16N_1988_tritium.min().min(),A16N_2003_tritium.min().min())
vmax = max(A16N_1988_tritium.max().max(),A16N_2003_tritium.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.875,wspace=0.05)
pm = axs[0,0].pcolormesh(A16N_1988_lat, A16N_1988_z, A16N_1988_tritium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='A16N tritium 1988 - observations')
pm = axs[0,1].pcolormesh(A16N_2003_lat, A16N_2003_z, A16N_2003_tritium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='A16N tritium 2003 - observations')
pm = axs[1,0].pcolormesh(model_A16N_1988.lat, -model_A16N_1988.z_t/100, model_tritium_1988, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='A16N tritium 1988 - model')
pm = axs[1,1].pcolormesh(model_A16N_2003.lat, -model_A16N_2003.z_t/100, model_tritium_2003, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='A16N tritium 2003 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('TU')

Helium

In [15]:
A16N_1988_helium, A16N_1988_lat, A16N_1988_z = get_woce_section('A16N',1988,'helium','latitude')
A16N_2003_helium, A16N_2003_lat, A16N_2003_z = get_woce_section('A16N',2003,'helium','latitude')

model_helium_1988 = model_A16N_1988.HELIUM4.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_helium_2003 = model_A16N_2003.HELIUM4.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_helium_1988[model_helium_1988<0] = 0
model_helium_2003[model_helium_2003<0] = 0

vmin = min(A16N_1988_helium.min().min(),A16N_2003_helium.min().min())
vmax = max(A16N_1988_helium.max().max(),A16N_2003_helium.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.875,wspace=0.05)
pm = axs[0,0].pcolormesh(A16N_1988_lat, A16N_1988_z, A16N_1988_helium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='A16N helium 1988 - observations')
pm = axs[0,1].pcolormesh(A16N_2003_lat, A16N_2003_z, A16N_2003_helium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='A16N helium 2003 - observations')
pm = axs[1,0].pcolormesh(model_A16N_1988.lat, -model_A16N_1988.z_t/100, model_helium_1988, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='A16N helium 1988 - model')
pm = axs[1,1].pcolormesh(model_A16N_2003.lat, -model_A16N_2003.z_t/100, model_helium_2003, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='A16N helium 2003 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Neon

In [16]:
A16N_1988_neon, A16N_1988_lat, A16N_1988_z = get_woce_section('A16N',1988,'neon','latitude')
A16N_2003_neon, A16N_2003_lat, A16N_2003_z = get_woce_section('A16N',2003,'neon','latitude')

model_neon_1988 = model_A16N_1988.NEON.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_neon_2003 = model_A16N_2003.NEON.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_neon_1988[model_neon_1988<0] = 0
model_neon_2003[model_neon_2003<0] = 0

vmin = min(A16N_1988_neon.min().min(),A16N_2003_neon.min().min())
vmax = max(A16N_1988_neon.max().max(),A16N_2003_neon.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.875,wspace=0.05)
pm = axs[0,0].pcolormesh(A16N_1988_lat, A16N_1988_z, A16N_1988_neon.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='A16N neon 1988 - observations')
pm = axs[0,1].pcolormesh(A16N_2003_lat, A16N_2003_z, A16N_2003_neon.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='A16N neon 2003 - observations')
pm = axs[1,0].pcolormesh(model_A16N_1988.lat, -model_A16N_1988.z_t/100, model_neon_1988, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='A16N neon 1988 - model')
pm = axs[1,1].pcolormesh(model_A16N_2003.lat, -model_A16N_2003.z_t/100, model_neon_2003, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='A16N neon 2003 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Helium 3

In [17]:
A16N_1988_delhe3, A16N_1988_lat, A16N_1988_z = get_woce_section('A16N',1988,'delhe3','latitude')
A16N_2003_delhe3, A16N_2003_lat, A16N_2003_z = get_woce_section('A16N',2003,'delhe3','latitude')
A16N_1988_helium3 = Ir * (1 + A16N_1988_delhe3/100) * A16N_1988_helium
A16N_2003_helium3 = Ir * (1 + A16N_2003_delhe3/100) * A16N_2003_helium

model_helium3_1988 = model_A16N_1988.HELIUM3.transpose().values / rho / 1.e+3 # pmol/m^3 -> nmol/kg
model_helium3_2003 = model_A16N_2003.HELIUM3.transpose().values / rho / 1.e+3 # pmol/m^3 -> nmol/kg
model_helium3_1988[model_helium3_1988<0] = 0
model_helium3_2003[model_helium3_2003<0] = 0

vmin = min(A16N_1988_helium3.min().min(),A16N_2003_helium3.min().min())
vmax = max(A16N_1988_helium3.max().max(),A16N_2003_helium3.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.85,wspace=0.05)
pm = axs[0,0].pcolormesh(A16N_1988_lat, A16N_1988_z, A16N_1988_helium3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='A16N $^3$He 1988 - observations')
pm = axs[0,1].pcolormesh(A16N_2003_lat, A16N_2003_z, A16N_2003_helium3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='A16N $^3$He 2003 - observations')
pm = axs[1,0].pcolormesh(model_A16N_1988.lat, -model_A16N_1988.z_t/100, model_helium3_1988, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='A16N $^3$He 1988 - model')
pm = axs[1,1].pcolormesh(model_A16N_2003.lat, -model_A16N_2003.z_t/100, model_helium3_2003, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='A16N $^3$He 2003 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Excess helium 3

In [18]:
A16N_1988_temp, _, _ = get_woce_section_ts('A16N',1988,'temperature','latitude')
A16N_1988_sal,  _, _ = get_woce_section_ts('A16N',1988,'salinity','latitude')
A16N_2003_temp, _, _ = get_woce_section_ts('A16N',2003,'temperature','latitude')
A16N_2003_sal,  _, _ = get_woce_section_ts('A16N',2003,'salinity','latitude')

A16N_1988_excess_he3 = excess_helium3(A16N_1988_temp, A16N_1988_sal, A16N_1988_helium3, None)
A16N_2003_excess_he3 = excess_helium3(A16N_2003_temp, A16N_2003_sal, A16N_2003_helium3, A16N_2003_neon)

model_excess_he3_1988 = (model_A16N_1988.HELIUM3 - model_A16N_1988.INERT_HELIUM3).transpose().values 
model_excess_he3_1988 = model_excess_he3_1988 / rho / 1.e+3 # pmol/m^3->nmol/kg
model_excess_he3_1988[model_excess_he3_1988<0] = 0
model_excess_he3_2003 = (model_A16N_2003.HELIUM3 - model_A16N_2003.INERT_HELIUM3).transpose().values 
model_excess_he3_2003 = model_excess_he3_2003 / rho / 1.e+3 # pmol/m^3->nmol/kg
model_excess_he3_2003[model_excess_he3_2003<0] = 0

vmin = min(A16N_1988_excess_he3.min().min(),A16N_2003_excess_he3.min().min())
vmax = max(A16N_1988_excess_he3.max().max(),A16N_2003_excess_he3.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.85,wspace=0.05)
pm = axs[0,0].pcolormesh(A16N_1988_lat, A16N_1988_z, A16N_1988_excess_he3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='A16N excess $^3$He 1988 - observations')
pm = axs[0,1].pcolormesh(A16N_2003_lat, A16N_2003_z, A16N_2003_excess_he3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='A16N excess $^3$He 2003 - observations')
pm = axs[1,0].pcolormesh(model_A16N_1988.lat, -model_A16N_1988.z_t/100, model_excess_he3_1988, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='A16N excess $^3$He 1988 - model')
pm = axs[1,1].pcolormesh(model_A16N_2003.lat, -model_A16N_2003.z_t/100, model_excess_he3_2003, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='A16N excess $^3$He 2003 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

P16N

Read model results

In [19]:
model_P16N_1991 = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/%s_P16N_1991.nc'%case)
model_P16N_2006 = xr.open_dataset('/home/ivan/Data/Postproc/Tritium-3H/%s_P16N_2006.nc'%case)

Tritium

In [20]:
P16N_1991_tritium, P16N_1991_lat, P16N_1991_z = get_woce_section('P16N',1991,'tritium','latitude')
P16N_2006_tritium, P16N_2006_lat, P16N_2006_z = get_woce_section('P16N',2006,'tritium','latitude')

model_tritium_1991 = model_P16N_1991.TRITIUM.transpose().values/0.10995 # pmol/m^3 -> TU
model_tritium_2006 = model_P16N_2006.TRITIUM.transpose().values/0.10995 # pmol/m^3 -> TU
model_tritium_1991[model_tritium_1991<0] = 0
model_tritium_2006[model_tritium_2006<0] = 0

vmin = min(P16N_1991_tritium.min().min(),P16N_2006_tritium.min().min())
vmax = max(P16N_1991_tritium.max().max(),P16N_2006_tritium.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.875,wspace=0.05)
pm = axs[0,0].pcolormesh(P16N_1991_lat, P16N_1991_z[5:], P16N_1991_tritium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='P16N tritium 1991 - observations')
pm = axs[0,1].pcolormesh(P16N_2006_lat, P16N_2006_z, P16N_2006_tritium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='P16N tritium 2006 - observations')
pm = axs[1,0].pcolormesh(model_P16N_1991.lat, -model_P16N_1991.z_t/100, model_tritium_1991, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='P16N tritium 1991 - model')
pm = axs[1,1].pcolormesh(model_P16N_2006.lat, -model_P16N_2006.z_t/100, model_tritium_2006, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='P16N tritium 2006 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('TU')

Helium

In [21]:
P16N_1991_helium, P16N_1991_lat, P16N_1991_z = get_woce_section('P16N',1991,'helium','latitude')
P16N_2006_helium, P16N_2006_lat, P16N_2006_z = get_woce_section('P16N',2006,'helium','latitude')

model_helium_1991 = model_P16N_1991.HELIUM4.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_helium_2006 = model_P16N_2006.HELIUM4.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_helium_1991[model_helium_1991<0] = 0
model_helium_2006[model_helium_2006<0] = 0

vmin = min(P16N_1991_helium.min().min(),P16N_2006_helium.min().min())
#vmax = max(P16N_1991_helium.max().max(),P16N_2006_helium.max().max())
vmax = P16N_2006_helium.max().max()

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.875,wspace=0.05)
pm = axs[0,0].pcolormesh(P16N_1991_lat, P16N_1991_z, P16N_1991_helium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='P16N helium 1991 - observations')
pm = axs[0,1].pcolormesh(P16N_2006_lat, P16N_2006_z, P16N_2006_helium.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='P16N helium 2006 - observations')
pm = axs[1,0].pcolormesh(model_P16N_1991.lat, -model_P16N_1991.z_t/100, model_helium_1991, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='P16N helium 1991 - model')
pm = axs[1,1].pcolormesh(model_P16N_2006.lat, -model_P16N_2006.z_t/100, model_helium_2006, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='P16N helium 2006 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Neon

In [22]:
P16N_1991_neon, P16N_1991_lat, P16N_1991_z = get_woce_section('P16N',1991,'neon','latitude')
P16N_2006_neon, P16N_2006_lat, P16N_2006_z = get_woce_section('P16N',2006,'neon','latitude')

model_neon_1991 = model_P16N_1991.NEON.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_neon_2006 = model_P16N_2006.NEON.transpose().values / rho * 1.e+3 # umol/m^3 -> nmol/kg
model_neon_1991[model_neon_1991<0] = 0
model_neon_2006[model_neon_2006<0] = 0

vmin = min(P16N_1991_neon.min().min(),P16N_2006_neon.min().min())
vmax = max(P16N_1991_neon.max().max(),P16N_2006_neon.max().max())

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.875,wspace=0.05)
pm = axs[0,0].pcolormesh(P16N_1991_lat, P16N_1991_z[9:], P16N_1991_neon.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='P16N neon 1991 - observations')
pm = axs[0,1].pcolormesh(P16N_2006_lat, P16N_2006_z, P16N_2006_neon.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='P16N neon 2006 - observations')
pm = axs[1,0].pcolormesh(model_P16N_1991.lat, -model_P16N_1991.z_t/100, model_neon_1991, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='P16N neon 1991 - model')
pm = axs[1,1].pcolormesh(model_P16N_2006.lat, -model_P16N_2006.z_t/100, model_neon_2006, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='P16N neon 2006 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Helium 3

In [23]:
P16N_1991_delhe3, P16N_1991_lat, P16N_1991_z = get_woce_section('P16N',1991,'delhe3','latitude')
P16N_2006_delhe3, P16N_2006_lat, P16N_2006_z = get_woce_section('P16N',2006,'delhe3','latitude')
P16N_1991_helium3 = Ir * (1 + P16N_1991_delhe3/100) * P16N_1991_helium
P16N_2006_helium3 = Ir * (1 + P16N_2006_delhe3/100) * P16N_2006_helium

model_helium3_1991 = model_P16N_1991.HELIUM3.transpose().values / rho / 1.e+3 # pmol/m^3 -> nmol/kg
model_helium3_2006 = model_P16N_2006.HELIUM3.transpose().values / rho / 1.e+3 # pmol/m^3 -> nmol/kg
model_helium3_1991[model_helium3_1991<0] = 0
model_helium3_2006[model_helium3_2006<0] = 0

vmin = min(P16N_1991_helium3.min().min(),P16N_2006_helium3.min().min())
#vmax = max(P16N_1991_helium3.max().max(),P16N_2006_helium3.max().max())
vmax = P16N_2006_helium3.max().max()

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.85,wspace=0.05)
pm = axs[0,0].pcolormesh(P16N_1991_lat, P16N_1991_z, P16N_1991_helium3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='P16N $^3$He 1991 - observations')
pm = axs[0,1].pcolormesh(P16N_2006_lat, P16N_2006_z, P16N_2006_helium3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='P16N $^3$He 2006 - observations')
pm = axs[1,0].pcolormesh(model_P16N_1991.lat, -model_P16N_1991.z_t/100, model_helium3_1991, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='P16N $^3$He 1991 - model')
pm = axs[1,1].pcolormesh(model_P16N_2006.lat, -model_P16N_2006.z_t/100, model_helium3_2006, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='P16N $^3$He 2006 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Excess helium 3

In [24]:
P16N_1991_temp, _, _ = get_woce_section_ts('P16N',1991,'temperature','latitude')
P16N_1991_sal,  _, _ = get_woce_section_ts('P16N',1991,'salinity','latitude')
P16N_2006_temp, _, _ = get_woce_section_ts('P16N',2006,'temperature','latitude')
P16N_2006_sal,  _, _ = get_woce_section_ts('P16N',2006,'salinity','latitude')

P16N_1991_excess_he3 = excess_helium3(P16N_1991_temp, P16N_1991_sal, P16N_1991_helium3, None)
P16N_2006_excess_he3 = excess_helium3(P16N_2006_temp, P16N_2006_sal, P16N_2006_helium3, P16N_2006_neon)

model_excess_he3_1991 = (model_P16N_1991.HELIUM3 - model_P16N_1991.INERT_HELIUM3).transpose().values 
model_excess_he3_1991 = model_excess_he3_1991 / rho / 1.e+3 # pmol/m^3->nmol/kg
model_excess_he3_1991[model_excess_he3_1991<0] = 0
model_excess_he3_2006 = (model_P16N_2006.HELIUM3 - model_P16N_2006.INERT_HELIUM3).transpose().values 
model_excess_he3_2006 = model_excess_he3_2006 / rho / 1.e+3 # pmol/m^3->nmol/kg
model_excess_he3_2006[model_excess_he3_2006<0] = 0

vmin = min(P16N_1991_excess_he3.min().min(),P16N_2006_excess_he3.min().min())
#vmax = max(P16N_1991_excess_he3.max().max(),P16N_2006_excess_he3.max().max())
vmax = P16N_2006_excess_he3.max().max()

fig, axs = plt.subplots(2,2,sharex=False,sharey=True,figsize=(9.5,9.6))
fig.subplots_adjust(left=0.1,right=0.85,wspace=0.05)
pm = axs[0,0].pcolormesh(P16N_1991_lat, P16N_1991_z, P16N_1991_excess_he3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,0].set(ylabel='depth (m)',title='P16N excess $^3$He 1991 - observations')
pm = axs[0,1].pcolormesh(P16N_2006_lat, P16N_2006_z, P16N_2006_excess_he3.values[:,1:], cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[0,1].set(title='P16N excess $^3$He 2006 - observations')
pm = axs[1,0].pcolormesh(model_P16N_1991.lat, -model_P16N_1991.z_t/100, model_excess_he3_1991, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,0].set(xlabel='latitude ($^{\circ}$N)',ylabel='depth (m)',title='P16N excess $^3$He 1991 - model')
pm = axs[1,1].pcolormesh(model_P16N_2006.lat, -model_P16N_2006.z_t/100, model_excess_he3_2006, cmap=plt.cm.plasma,
                      vmin=vmin, vmax=vmax)
_ = axs[1,1].set(xlabel='latitude ($^{\circ}$N)',title='P16N excess $^3$He 2006 - model')
_ = axs[1,1].set_ylim(zbins.min(), zbins.max())
l, b, w, h = axs[0,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('nmol kg$^{-1}$')

Save WOCE section positions to file

In [25]:
def get_section_location(sec,year):
    data = database[(database.sect_id==sec) & (database.stadate.dt.year==year)][info_cols]
    data = data.drop_duplicates('latitude').sort_values('latitude')
    return np.c_[data.longitude.values, data.latitude.values]

A16N_1988_lonlat = get_section_location('A16N', 1988)
A16N_2003_lonlat = get_section_location('A16N', 2003)
P16N_1991_lonlat = get_section_location('P16N', 1991)
P16N_2006_lonlat = get_section_location('P16N', 2006)

np.savez('../data/woce_sec_lonlat.npz',
         A16N_1988_lonlat = A16N_1988_lonlat,
         A16N_2003_lonlat = A16N_2003_lonlat,
         P16N_1991_lonlat = P16N_1991_lonlat,
         P16N_2006_lonlat = P16N_2006_lonlat)