Compute $^3\!$He air-sea fluxes using the Stanley et al 2009, 2015 model¶
Created by Ivan Lima on Fri, 03 Jun 2016 11:15:18 -0400
Compute air-sea $^3\!$He fluxes using Stanley 2009, 2015 model, $u_{10}$ and sea level pressure from NCEP and surface $^3\!$He concentrations from a CESM test run initialized with $^3\!$He concentrations at saturation.
In [45]:
import xarray as xr
import pandas as pd
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import glob, os, datetime, netCDF4, gsw, warnings
from cartopy.util import add_cyclic_point
from scipy.interpolate import griddata
from mpl_utils import spd, center_cmap
from tr3he import *
%matplotlib inline
In [46]:
def fixdates(ds):
""" Convert netcdftime to datetime dates. ds = xarray DataSet """
units = ds['time'].units.replace('0000','1800')
dates = netCDF4.num2date(ds['time'] - 1,units=units,calendar=ds['time'].calendar)
dates = [datetime.datetime(*d.timetuple()[:3]) for d in dates]
ds['time'] = pd.DatetimeIndex(dates)
ds['time'].attrs['calendar'] = 'noleap'
def remap(da,vname,lname,units):
""" Interpolate NCEP DataArray (da) to CESM grid """
lon, lat, data = da.lon.values, da.lat.values, da.values
data, lon = add_cyclic_point(data,coord=lon) # add cyclic points to data & lon to avoid gaps
lon, lat = np.meshgrid(lon, lat)
newdata = griddata((lon.ravel(),lat.ravel()),data.ravel(),(dsin.TLONG,dsin.TLAT),method='linear')
newda = xr.DataArray(newdata.reshape((1,)+newdata.shape),coords=dsin.HMXL.coords,dims=dsin.HMXL.dims,name=vname)
newda = newda.where(dsin.REGION_MASK>0) # mask land & marginal seas
newda.attrs['long_name'] = lname
newda.attrs['units'] = units
return newda
def da2ma(da):
"""Convert xr.DataArray to numpy masked array"""
return np.ma.masked_where(np.isnan(da.data),da.data)
proj = ccrs.Robinson(central_longitude=210)
warnings.filterwarnings('ignore')
Read input data from model test run¶
In [47]:
datadir = '/bali/data/ilima/cesm/archive/test_GIAF_tr3he.gx1.004/ocn/hist/'
filelist = sorted(glob.glob(datadir + 'test_GIAF_tr3he.gx1.004.pop.h.????-??.nc'))
dsin = xr.open_dataset(filelist[0],decode_times=False)
fixdates(dsin)
dsin['PD'].values = dsin.PD.where(dsin.REGION_MASK>0) * 1.e+3 # mask marginal seas & g/cm^3 -> km/m^3
dsin['PD'].attrs['units'] = 'kg/m^3'
dsin['TEMP'].values = dsin.TEMP.where(dsin.REGION_MASK>0) # mask marginal seas
dsin['SALT'].values = dsin.SALT.where(dsin.REGION_MASK>0) # mask marginal seas
dsin['HELIUM3'].values = dsin.HELIUM3.where(dsin.REGION_MASK>0) # mask marginal seas
Get wind speed ($u_{10}$) and sea level pressure from NCEP¶
In [48]:
duvel = xr.open_dataset('/bali/data/ilima/cesm/cesm1_input/ocn/iaf/ncep.u_10.T62.1948.nc')
dvvel = xr.open_dataset('/bali/data/ilima/cesm/cesm1_input/ocn/iaf/ncep.v_10.T62.1948.nc')
dslp = xr.open_dataset('/bali/data/ilima/cesm/cesm1_input/ocn/iaf/ncep.slp_.T62.1948.nc')
u_ncep, v_ncep = duvel.isel(time=0), dvvel.isel(time=0)
slp_ncep = dslp.isel(time=0)
u10_ncep = np.sqrt(u_ncep.u_10**2 + v_ncep.v_10**2) # compute u10 from uvel & vvel
Remap wind speed ($u_{10}$) and sea level pressure to CESM grid¶
In [49]:
u10 = remap(u10_ncep,'u10','wind velocity at 10m','m/s')
slp = remap(slp_ncep.slp_,'slp','sea level pressure','Pa')
ds = xr.Dataset({'slp':slp, 'u10':u10})
ds = ds.merge(dsin[['REGION_MASK','HMXL','PD','TEMP','SALT','HELIUM3']])
ds.HELIUM3.values = ds.HELIUM3 * 1.e-12 # pmol/m^3 -> mol/m^3
ds.HELIUM3.attrs['units'] = 'mol/m^3'
Input fields maps¶
In [50]:
fig, axs = plt.subplots(nrows=2,ncols=2,subplot_kw={'projection':proj},figsize=(12,8))
fig.subplots_adjust(wspace=0.05,hspace=0.15)
axs = axs.ravel()
varlist = ['TEMP','SALT','u10','slp']
for n in range(len(varlist)):
vname = varlist[n]
lat, lon = ds.TLAT.values, ds.TLONG.values
latc, lonc = np.hstack((lat,lat[:,0,np.newaxis])), np.hstack((lon,lon[:,0,np.newaxis]))
if ds[vname].ndim == 4:
data = ds[vname][0,0,...]
else:
data = ds[vname][0,...]
data = da2ma(data)
axs[n].set_global()
axs[n].add_feature(cfeature.LAND,facecolor='#b0b0b0')
axs[n].set_title('%s (%s)'%(ds[vname].long_name,vname))
pm = axs[n].pcolormesh(lonc,latc,data,transform=ccrs.PlateCarree())
cb = fig.colorbar(pm,ax=axs[n],orientation='horizontal',pad=0.03)
cb.set_label(ds[vname].units)
if vname == 'slp':
cb.ax.set_xticklabels(cb.ax.get_xticklabels(),rotation=20,horizontalalignment='right')
$^3\!$He equilibrium concentration at standard atmospheric pressure¶
The $^3\!$He equilibrium concentration $\left(C_{eq}\right)$ at standard atmospheric pressure is given by:
$$ C_{eq} = S\,P_{atm}\,X_{He}\,I_r\,\alpha_S $$$S$ = $^4\!$He solubility (mol m$^{-3}$ Pa$^{-1}$), $P_{atm}$ = standard atmospheric pressure (101325 Pa), $X_{He}$ = $^4\!$He atmospheric mole fraction, $I_r$ = $^3\!$He/$^4\!$He isotopic ratio and $\alpha_S$ = temperature-dependent solubility fractionation factor.
In [51]:
sst, sss = ds.TEMP.isel(z_t=0), ds.SALT.isel(z_t=0)
he4_sol = comp_henry_he_sol_0(sst,sss)
alpha_sol = alpha_sol_he(sst)
he4_ceq = he4_sol * Xhe * Patm # 4He equilibrium concentration
he3_ceq = he4_ceq * Ir * alpha_sol # 3He equilibrium concentration
he4_ceq.attrs = {'long_name':'helium 4 equilibriun concentration','units':'mol/m^3'}
he3_ceq.attrs = {'long_name':'helium 3 equilibriun concentration','units':'mol/m^3'}
ds = ds.merge({'he4_ceq':he4_ceq, 'he3_ceq':he3_ceq})
Maps of surface $^4\!$He and $^3\!$He saturation concentrations¶
In [52]:
fig, axs = plt.subplots(ncols=2,subplot_kw={'projection':proj},figsize=(12,5))
fig.subplots_adjust(wspace=0.05)
axs = axs.ravel()
for ax in axs:
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
data = da2ma(he4_ceq.isel(time=0))
data = np.ma.masked_where(np.isnan(data),data) * 1.e+6 # mol/m^3 -> umol/m^3
pm = axs[0].pcolormesh(lonc,latc,data,transform=ccrs.PlateCarree())
axs[0].set_title(r'$^4$He equilibrium concentration')
cb = fig.colorbar(pm,ax=axs[0],orientation='horizontal',pad=0.03)
cb.set_label(r'$\mu$mol m$^{-3}$')
data = da2ma(he3_ceq.isel(time=0))
data = np.ma.masked_where(np.isnan(data),data) * 1.e+12 # mol/m^3 -> pmol/m^3
pm = axs[1].pcolormesh(lonc,latc,data,transform=ccrs.PlateCarree())
axs[1].set_title(r'$^3$He equilibrium concentration ($C_{eq}$)')
cb = fig.colorbar(pm,ax=axs[1],orientation='horizontal',pad=0.03)
cb.set_label(r'pmol m$^{-3}$')
Maps of surface $^3\!$He concentration and saturation¶
In [53]:
loncp = np.where(np.greater_equal(lon,lon[:,0].min()),lon-360,lon) # monotonic lon for contour plots
fig, axs = plt.subplots(ncols=2,subplot_kw={'projection':proj},figsize=(12,5))
fig.subplots_adjust(wspace=0.05)
axs = axs.ravel()
for ax in axs:
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
he3 = da2ma(ds.HELIUM3.isel(time=0,z_t=0)) * 1.e+12 # mol/m^3 -> pmol/m^3
pm = axs[0].pcolormesh(lonc,latc,he3,transform=ccrs.PlateCarree())
axs[0].set_title(r'surface $^3$He concentration ($C$)')
cb = fig.colorbar(pm,ax=axs[0],orientation='horizontal',pad=0.03)
cb.set_label(r'$\mu$mol m$^{-3}$')
eq = da2ma(he3_ceq.isel(time=0)) * 1.e+12 # mol/m^3 -> pmol/m^3
sat = (he3 - eq) / eq * 100
pm = axs[1].pcolormesh(lonc,latc,sat,transform=ccrs.PlateCarree(),cmap=plt.cm.seismic)
center_cmap(pm)
cl = axs[1].contour(loncp,ds.TLAT,sat,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
axs[1].set_title(r'$\frac{C-C_{eq}}{C_{eq}} \times 100$')
cb = fig.colorbar(pm,ax=axs[1],orientation='horizontal',pad=0.03)
cb.set_label(r'\%')
$^3\!$He gas fluxes¶
Diffusive gas exchange¶
The difusive gas exchange in mol m$^{-2}$ s$^{-1}$ is given by:
$$ F_D = K_s \left(C_{eq} - C\right) $$$K_s$ is the gas transfer velocity (m s$^{-1}$):
$$ K_s = 0.7 \times 8.6 \times 10^{-7}\;u_{10}^2 \sqrt{\frac{666}{\frac{1}{\alpha_D}S_C}} $$$u_{10}$ = wind speed at 10 m, $C$ = $^3\!$He concentration in the water. The Schmidt number ($S_c$) is scaled by the $^3$He/$^4$He diffusivity ratio ($\alpha_D$) from Bourg & Sposito (2008).
In [54]:
u10 = ds.u10
C = ds.HELIUM3.isel(z_t=0)
Ceq = ds.he3_ceq
Sc = comp_he_schmidt(sst) / alpha_diff # scale 4He Schmidt number by 3He/4He diffusivity ratio
Ks = 0.7 * 8.6e-7 * u10**2 * np.sqrt(660/Sc) # compute gas transfer velocity in m/s
Fd = Ks * (Ceq - C) # compute gas flux in mol/m^2/s
Fd.attrs = {'long_name':'diffusive gas exchange','units':'mol/m^2/s'}
ds = ds.merge({'Fd':Fd})
Plot map
In [55]:
fig, ax = plt.subplots(subplot_kw={'projection':proj},figsize=(8,6))
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
Fd_pmol = da2ma(Fd.isel(time=0)) * spd * 1.e+12 # mol/m^2/s -> pmol/m^2/d
ax.set_title(r'Diffusive gas exchange ($F_D$)')
pm = ax.pcolormesh(lonc,latc,Fd_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.seismic)
center_cmap(pm)
cl = ax.contour(loncp,ds.TLAT,Fd_pmol,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
cb = fig.colorbar(pm,ax=ax,orientation='horizontal',pad=0.03)
cb.set_label(r'pmol m$^{-2}$ d$^{-1}$')
vmin, vmax = pm.get_clim()
print('min = %.3g, max = %.3g'%(vmin, vmax))
Flux due to partially dissolved bubbles¶
The $^3\!$He flux due to partially trapped bubbles (mol m$^{-2}$ s$^{-1}$) is computed as:
$$F_P = 2.2 \times 10^{-3} \left(u_{10} - 2.27\right)^3 \beta\,D^{\frac{2}{3}} \frac{P_b - P_w}{R\,T}$$$\beta$ is the Bunsen solubility coefficient for $^3\!$He:
$$ \beta = \alpha_S\,\beta_{^4\!He} $$$\beta_{^4\!He}$ = Bunsen solubility for $^4\!$He, $\alpha_S$ = temperature dependent solubility fractionation factor.
$D$ is the diffudivity coefficient for $^3\!$He (m$^2$ s$^{-1}$):
$$ D = \alpha_D\,D_{^4\!He} $$$D_{^4\!He}$ = diffudivity coefficient for $^4\!$He, $\alpha_D$ = $^3\!$He/$^4\!$He diffusivity ratio.
$P_b$ is $^3\!$He partial pressure in the bubble:
$$ P_b = X_{He} I_r \left(P_{slp} + \rho g z_b\right) $$$X_{He}$ = $^4\!$He atmospheric mole fraction, $I_r$ = $^3\!$He/$^4\!$He isotopic ratio, $P_{slp}$ = sea level pressure (Pa), $\rho$ = density of seawater (kg m$^{-3}$), $g$ = gravitational acceleration (m s$^{-2}$) and $z_b$ = depth of the bubble (m). The bubble depth is a function of $u_{10}$ and computed as:
$$ z_b = 0.15 u_{10} - 0.55 $$$P_w$ is the $^3\!$He partial pressure in the water:
$$ P_w = \frac{C}{\alpha_S S} $$$C$ = $^3\!$He concentration in the water, $S$ = $^4\!$He solubility (mol m$^{-3}$ Pa$^{-1}$), and $\alpha_S$ = temperature dependent solubility fractionation factor.
In [56]:
rhow = ds.PD.isel(z_t=0)
Pslp = ds.slp
C = ds.HELIUM3.isel(z_t=0)
B = comp_bunsen_he_sol_0(sst,sss) * alpha_sol # Bunsen solubility coefficient
D = comp_he_diff(sst) * alpha_diff # diffusivity coefficient
zb = 0.15 * u10 - 0.55 # bubble depth
zb.data = np.where(zb<0,0,zb)
Pb = Xhe * Ir * (Pslp + rhow * g * zb) # 3He partial pressure in bubble
Pw = C / (he4_sol * alpha_sol) # 3He partial pressure in water
uwork = np.where(u10>2.27,u10-2.27,0)
Fp = 2.2e-3 * (uwork)**3 * B * D**(2/3) * (Pb - Pw)/(R*(sst+T0K)) # mol/m^2/s
Fp.attrs = {'long_name':'flux due to partially dissolved bubbles','units':'mol/m^2/s'}
ds = ds.merge({'Fp':Fp})
Plot map
In [57]:
fig, ax = plt.subplots(subplot_kw={'projection':proj},figsize=(8,6))
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
Fp_pmol = da2ma(Fp.isel(time=0)) * spd * 1.e+12 # mol/m^2/s -> pmol/m^2/d
ax.set_title(r'flux due to partially dissolved bubbles ($F_P$)')
pm = ax.pcolormesh(lonc,latc,Fp_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.seismic)
center_cmap(pm)
#cl = ax.contour(loncp,ds.TLAT,Fp_pmol,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
cb = fig.colorbar(pm,ax=ax,orientation='horizontal',pad=0.03)
cb.set_label(r'pmol m$^{-2}$ d$^{-1}$')
vmin, vmax = pm.get_clim()
print('min = %.3g, max = %.3g'%(vmin, vmax))
Flux due to completely dissolved bubbles¶
The $^3\!$He flux due to completely trapped bubbles (mol m$^{-2}$ s$^{-1}$) is given by:
$$ F_C = 9.1 \times 10^{-11} \left(u_{10} - 2.27\right)^3 \frac{P_{a}}{R\,T} $$$P_{a}$ is the partial pressure of $^3\!$He in the atmosphere:
$$ P_{a} = X_{He} P_{slp} I_r $$$R$ = ideal gas constant and $T$ = temperature in Kelvin, $X_{He}$ = $^4$He atmospheric mole fraction, $P_{slp}$ = sea level pressure (Pa) and $I_r$ = $^3$He/$^4$He isotopic ratio.
In [58]:
Pa = Xhe * Pslp * Ir
Fc = 9.1e-11 * (uwork)**3 * Pa/(R*(sst+T0K))
Fc.attrs = {'long_name':'flux due to completely dissolved bubbles','units':'mol/m^2/s'}
ds = ds.merge({'Fc':Fc})
Plot map
In [59]:
fig, ax = plt.subplots(subplot_kw={'projection':proj},figsize=(8,6))
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
Fc_pmol = da2ma(Fc.isel(time=0)) * spd * 1.e+12 # mol/m^2/s -> pmol/m^2/d
ax.set_title(r'flux due to completely dissolved bubbles ($F_C$)')
pm = ax.pcolormesh(lonc,latc,Fc_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.seismic)
center_cmap(pm)
cb = fig.colorbar(pm,ax=ax,orientation='horizontal',pad=0.03)
#cl = ax.contour(loncp,ds.TLAT,Fc_pmol,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
cb.set_label(r'pmol m$^{-2}$ d$^{-1}$')
vmin, vmax = pm.get_clim()
print('min = %.3g, max = %.3g'%(vmin, vmax))
Total gas flux¶
In [60]:
Ft = Fd + Fp + Fc
fig, ax = plt.subplots(subplot_kw={'projection':proj},figsize=(8,6))
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
Ft_pmol = da2ma(Ft.isel(time=0)) * spd * 1.e+12 # mol/m^2/s -> pmol/m^2/d
ax.set_title(r'total gas flux ($F_{tot}$)')
pm = ax.pcolormesh(lonc,latc,Ft_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.seismic)
center_cmap(pm)
cb = fig.colorbar(pm,ax=ax,orientation='horizontal',pad=0.03)
cl = ax.contour(loncp,ds.TLAT,Ft_pmol,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
cb.set_label(r'pmol m$^{-2}$ d$^{-1}$')
vmin, vmax = pm.get_clim()
print('min = %.3g, max = %.3g'%(vmin, vmax))
