Diagnostic plots: time series

Created by Ivan Lima on Thu Jul 20 2017 15:47:57 -0400

In [1]:
%matplotlib notebook
import pandas as pd
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import os
from matplotlib import colors
from cesm_utils import *
from cases import *
from datetime import datetime
_ = plt.rcParams['xtick.minor.bottom'] = False # fix strange minor ticks in time series plots
print('Last updated on %s'%datetime.now().ctime())

Last updated on Tue Apr 30 09:15:47 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 data

In [3]:
datadir = '/home/ivan/Data/Postproc/Tritium-3H'
store = pd.HDFStore(os.path.join(datadir,'%s_tseries.h5'%case))
df_tseries = store['df']
store.close()
df_annual = df_tseries.groupby(df_tseries.index.year).mean() # compute annual means

$^3\!$H inventory & fluxes

Inventory & fluxes

In [4]:
fig, axs = plt.subplots(nrows=2,ncols=2,sharex=True,figsize=(9.5,7.125))
fig.subplots_adjust(wspace=0.25)
fig.delaxes(axs[1,1])
_ = df_tseries.TRITIUM.plot(ax=axs[0,0],title='$^3\!$H inventory')
_ = df_tseries.TRITIUM_DECAY.plot(ax=axs[0,1],title='$^3\!$H decay')
_ = df_tseries.STF_TRITIUM.plot(ax=axs[1,0],title='$^3\!$H total surface flux')
_ = axs[0,0].set_ylabel(r'mol')
_ = axs[0,1].set_ylabel(r'mol d$^{-1}$')
_ = axs[1,0].set_ylabel(r'mol d$^{-1}$')

Surface fluxes

In [5]:
fig, axs = plt.subplots(nrows=3,ncols=2,sharex=True,figsize=(9.5,9))
fig.delaxes(axs[2,1])
fig.subplots_adjust(wspace=0.25)
_ = df_tseries.STF_TRITIUM_PREC.plot(ax=axs[0,0],title='$^3\!$H precipitation flux')
_ = df_tseries.STF_TRITIUM_EVAP.plot(ax=axs[0,1],title='$^3\!$H evaporation flux')
_ = df_tseries.STF_TRITIUM_EVAP_UP.plot(ax=axs[1,0],title='upward $^3\!$H vapor flux')
_ = df_tseries.STF_TRITIUM_EVAP_DOWN.plot(ax=axs[1,1],title='downward $^3\!$H vapor flux')
_ = df_tseries.STF_TRITIUM_ROFF.plot(ax=axs[2,0],title='$^3\!$H river flux')
for ax in axs.ravel():
    _ = ax.set_ylabel(r'mol d$^{-1}$')

Precipitation & evaporation fluxes

In [6]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,3.4))
fig.subplots_adjust(wspace=0.3)
_ = df_tseries.PREC_F.plot(ax=ax1,title='precipitation')
_ = df_tseries.EVAP_F.plot(ax=ax2,title='evaporation')
_ = ax1.set_ylabel(r'm$^{-3}$ d$^{-1}$')
_ = ax2.set_ylabel(r'm$^{-3}$ d$^{-1}$')

Compare changes in $^3\!$H inventory with surface flux minus decay

In [7]:
dpm2 = np.resize(dpm,len(df_tseries)) # days per month
dtritium_dt = df_tseries.TRITIUM.diff()/dpm2     # mols/day
df_tseries['SMS_TRITIUM'] = df_tseries.STF_TRITIUM - df_tseries.TRITIUM_DECAY
fig, axs = plt.subplots(nrows=2,ncols=2,sharex=True,figsize=(9.5,7.125))
fig.subplots_adjust(wspace=0.25)
_ = dtritium_dt.plot(ax=axs[0,0],title=r'$\frac{\Delta ^3\!H}{\Delta t}$')
_ = df_tseries.STF_TRITIUM.plot(ax=axs[0,1],title='$^3\!H$ surface flux')
_ = df_tseries.TRITIUM_DECAY.plot(ax=axs[1,0],title='$^3\!H$ decay')
_ = df_tseries.SMS_TRITIUM.plot(ax=axs[1,1],title='$^3\!H$ surface flux - $^3\!H$ decay')
for ax in axs.ravel():
    _ = ax.set_ylabel(r'mol d$^{-1}$')
    
# df = pd.DataFrame({'$\Delta$3H/$\Delta$t':dtritium_dt,'flux - decay':df_tseries.SMS_TRITIUM})
# df.tail(12)

Model's global mean $^3\!$H flux

In [8]:
mean_flux = df_tseries.STF_TRITIUM.mean()
print('mean tritium flux = %.6f mol/day'%mean_flux)

mean tritium flux = 6.091047 mol/day

Estimated $^3\!$H decay based on inventory in last month of simulation

In [9]:
hlife = 12.31 * 365 # days
lam = -np.log(0.5)/hlife
decay = lam * df_tseries.TRITIUM[-1]
print('decay = %.6f mol/day'%decay)

decay = 1.746068 mol/day

Model's $^3\!$H decay in last month of simulation

In [10]:
print('decay = %.6f mol/day'%df_tseries.TRITIUM_DECAY[-1])

decay = 1.891358 mol/day

$^3\!$He inventory & air-sea fluxes

CTB = completely trapped bubbles, PTB = partially trapped bubbles, DGE = difusive gas exchange

In [11]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,4))
helium = df_tseries.HELIUM3 * 1.e-6 # mol -> 10^6 mol
_ = helium.plot(ax=ax1,title=r'$^3\!$He inventory')
_ = ax1.set_ylabel(r'$10^6$ mol')
_ = df_tseries[['STF_HELIUM3','STF_HELIUM3_CTB','STF_HELIUM3_PTB',
                'STF_HELIUM3_DGE']].plot(ax=ax2,title=r'$^3\!$He air-sea flux')
_ = ax2.hlines(0,df_tseries.index[0],df_tseries.index[-1],linestyles='dashed',linewidth=1)
_ = ax2.set_ylabel(r'mol d$^{-1}$')

Global mean fluxes (mol d$^{-1}$)

In [12]:
df_tseries[['STF_HELIUM3','STF_HELIUM3_CTB','STF_HELIUM3_PTB','STF_HELIUM3_DGE']].mean()
Out[12]:
STF_HELIUM3        -4.499278
STF_HELIUM3_CTB    50.240875
STF_HELIUM3_PTB    -1.065948
STF_HELIUM3_DGE   -53.674205
dtype: float64

Annual means

CTB = completely trapped bubbles, PTB = partially trapped bubbles, DGE = difusive gas exchange

In [13]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,4))
helium_annual = df_annual['HELIUM3'] * 1.e-6 # mol -> 10^6 mol
_ = helium_annual.plot(ax=ax1,title=r'$^3\!$He inventory')
_ = ax1.set_ylabel(r'$10^6$ mol')
_ = df_annual[['STF_HELIUM3','STF_HELIUM3_CTB','STF_HELIUM3_PTB',
               'STF_HELIUM3_DGE']].plot(ax=ax2,title=r'$^3\!$He air-sea flux')
_ = ax2.hlines(0,df_annual.index[0],df_annual.index[-1],linestyles='dashed',linewidth=1)
_ = ax2.set_ylabel(r'mol d$^{-1}$')

Inert $^3\!$He inventory & air-sea fluxes

CTB = completely trapped bubbles, PTB = partially trapped bubbles, DGE = difusive gas exchange

In [14]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,4))
helium = df_tseries[['INERT_HELIUM3']] * 1.e-6 # mol -> 10^6 mol
_ = helium.plot(ax=ax1,title=r'inert $^3\!$He inventory')
_ = ax1.set_ylabel(r'$10^6$ mol')
_ = df_tseries[['STF_INERT_HELIUM3','STF_INERT_HELIUM3_CTB','STF_INERT_HELIUM3_PTB',
                'STF_INERT_HELIUM3_DGE']].plot(ax=ax2,title=r'inert $^3\!$He air-sea flux')
_ = ax2.hlines(0,df_tseries.index[0],df_tseries.index[-1],linestyles='dashed',linewidth=1)
_ = ax2.set_ylabel(r'mol d$^{-1}$')

Global mean fluxes (mol d$^{-1}$)

In [15]:
df_tseries[['STF_INERT_HELIUM3','STF_INERT_HELIUM3_CTB','STF_INERT_HELIUM3_PTB','STF_INERT_HELIUM3_DGE']].mean()
Out[15]:
STF_INERT_HELIUM3        -0.129006
STF_INERT_HELIUM3_CTB    50.240875
STF_INERT_HELIUM3_PTB     1.016437
STF_INERT_HELIUM3_DGE   -51.386317
dtype: float64

Compare changes in $^3\!$He inventory with air-sea flux plus input from $^3\!$H decay

Data is very noisy, so we're looking at cumulative changes and fluxes.

In [16]:
cum_dhelium3 = (df_tseries.HELIUM3.diff()).cumsum() # mol
cum_flux_plus_decay = ((df_tseries.STF_HELIUM3 + df_tseries.TRITIUM_DECAY) * dpm2).cumsum() # mol
df = pd.DataFrame({r'$\sum$ $\Delta$3He':cum_dhelium3,r'$\sum$ (flux + 3H decay)$\times \Delta$t':cum_flux_plus_decay})
fig, ax = plt.subplots()
_ = df.plot(ax=ax)
_ = ax.set_ylabel(r'mol')
df.tail()
Out[16]:
$\sum$ $\Delta$3He $\sum$ (flux + 3H decay)$\times \Delta$t
2009-08-31 27821.841440 27888.718970
2009-09-30 27867.659537 27969.887756
2009-10-31 27857.739393 27863.035956
2009-11-30 27827.636732 27812.742461
2009-12-31 27679.769425 27663.223916

Excess $^3\!$He

In [17]:
df_tseries['EXCESS_HELIUM3'] = (df_tseries.HELIUM3 - df_tseries.INERT_HELIUM3)
df_tseries['STF_EXCESS_HELIUM3'] = df_tseries.STF_HELIUM3 - df_tseries.STF_INERT_HELIUM3
fig, axs = plt.subplots(ncols=2,nrows=2,sharex=True,figsize=(9.5,7.125))
axs = axs.ravel()
fig.delaxes(axs[-1])
fig.subplots_adjust(wspace=0.25)
_ = df_tseries.EXCESS_HELIUM3.plot(ax=axs[0],title=r'excess $^3\!$He inventory')
_ = axs[0].set_ylabel(r'mol')
_ = df_tseries.STF_EXCESS_HELIUM3.plot(ax=axs[1],title=r'excess $^3\!$He air-sea flux')
_ = axs[1].set_ylabel(r'mol d$^{-1}$')
_ = df_tseries.TRITIUM_DECAY.plot(ax=axs[2],title=r'$^3\!$H decay')
_ = axs[2].set_ylabel(r'mol d$^{-1}$')

$^4\!$He inventory & air-sea fluxes

CTB = completely trapped bubbles, PTB = partially trapped bubbles, DGE = difusive gas exchange

In [18]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,4))
helium4 = df_tseries.HELIUM4 * 1.e-12 # mol -> 10^12 mol
_ = helium4.plot(ax=ax1,title=r'$^4\!$He inventory')
_ = ax1.set_ylabel(r'$10^{12}$ mol')
_ = df_tseries[['STF_HELIUM4','STF_HELIUM4_CTB','STF_HELIUM4_PTB',
                'STF_HELIUM4_DGE']].plot(ax=ax2,title=r'$^4\!$He air-sea flux')
_ = ax2.hlines(0,df_tseries.index[0],df_tseries.index[-1],linestyles='dashed',linewidth=1)
_ = ax2.set_ylabel(r'mol d$^{-1}$')

Global mean fluxes (mol d$^{-1}$)

In [19]:
df_tseries[['STF_HELIUM4','STF_HELIUM4_CTB','STF_HELIUM4_PTB','STF_HELIUM4_DGE']].mean()
Out[19]:
STF_HELIUM4       -1.354277e+05
STF_HELIUM4_CTB    3.630121e+07
STF_HELIUM4_PTB    6.080421e+05
STF_HELIUM4_DGE   -3.704468e+07
dtype: float64

Neon inventory & air-sea fluxes

CTB = completely trapped bubbles, PTB = partially trapped bubbles, DGE = difusive gas exchange

In [20]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,4))
fig.subplots_adjust(wspace=0.25)
neon = df_tseries.NEON * 1.e-12 # mol -> 10^12 mol
_ = neon.plot(ax=ax1,title=r'Ne inventory')
_ = ax1.set_ylabel(r'$10^{12}$ mol')
_ = df_tseries[['STF_NEON','STF_NEON_CTB','STF_NEON_PTB','STF_NEON_DGE']].plot(ax=ax2,title=r'Ne air-sea flux')
_ = ax2.hlines(0,df_tseries.index[0],df_tseries.index[-1],linestyles='dashed',linewidth=1)
_ = ax2.set_ylabel(r'mol d$^{-1}$')

Global mean fluxes (mol d$^{-1}$)

In [21]:
df_tseries[['STF_NEON','STF_NEON_CTB','STF_NEON_PTB','STF_NEON_DGE']].mean()
Out[21]:
STF_NEON       -8.447269e+05
STF_NEON_CTB    1.259458e+08
STF_NEON_PTB   -4.206082e+06
STF_NEON_DGE   -1.225844e+08
dtype: float64

Argon inventory & air-sea fluxes

CTB = completely trapped bubbles, PTB = partially trapped bubbles, DGE = difusive gas exchange

In [22]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,4))
fig.subplots_adjust(wspace=0.25)
argon = df_tseries.ARGON * 1.e-15 # mol -> 10^15 mol
_ = argon.plot(ax=ax1,title=r'Ar inventory')
_ = ax1.set_ylabel(r'$10^{15}$ mol')
_ = df_tseries[['STF_ARGON','STF_ARGON_CTB','STF_ARGON_PTB','STF_ARGON_DGE']].plot(ax=ax2,title=r'Ne air-sea flux')
_ = ax2.hlines(0,df_tseries.index[0],df_tseries.index[-1],linestyles='dashed',linewidth=1)
_ = ax2.set_ylabel(r'mol d$^{-1}$')

Global mean fluxes (mol d$^{-1}$)

In [23]:
df_tseries[['STF_ARGON','STF_ARGON_CTB','STF_ARGON_PTB','STF_ARGON_DGE']].mean()
Out[23]:
STF_ARGON       -2.591903e+09
STF_ARGON_CTB    6.470483e+10
STF_ARGON_PTB    4.444829e+10
STF_ARGON_DGE   -1.117450e+11
dtype: float64

$^3\!$H concentration in precipitation (GNIP)

In [24]:
infile = os.path.join(datadir,'%s_tritium_prec.nc'%case)
ds = xr.open_dataset(infile,decode_times=False)
ds = fixtime(ds,1824)
ds['time'] = ds.time.to_index().to_datetimeindex()
In [25]:
tlat = ds.TLAT[:-13]
data = da2ma(ds.TRITIUM_PREC[:,:-13])
xt = pd.date_range('1950-01-01','2010-12-31',freq='10A')
fig, ax = plt.subplots()
pm = ax.pcolormesh(ds.time,tlat,data.transpose(),norm=colors.LogNorm(),cmap=plt.cm.YlOrRd)
_ = ax.set(xticks=xt,xlabel='time',ylabel='latitude',title=r'$^3\!$H concentration in precipitation')
cb = fig.colorbar(pm,ax=ax)
cb.set_label(r'pmol m$^{-3}$')

$^3\!$H concentration in rivers (GNIP)

In [26]:
data = da2ma(ds.TRITIUM_RIVER[:,:-13])
fig, ax = plt.subplots()
pm = ax.pcolormesh(ds.time,tlat,data.transpose(),norm=colors.LogNorm(),cmap=plt.cm.YlOrRd)
_ = ax.set(xticks=xt,xlabel='time',ylabel='latitude',title=r'$^3\!$H concentration in rivers')
cb = fig.colorbar(pm,ax=ax)
cb.set_label(r'pmol m$^{-3}$')