Compute $^3\!$He air-sea fluxes using the Liang et al 2013 model

Created by Ivan Lima on Wed, 14 Sep 2016 09:38:26 -0400

Compute air-sea $^3\!$He fluxes using Liang et al 2013 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 [1]:
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 notebook
In [2]:
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 calc_cd10(u10):
    """ Compute drag coefficient from u10 (m/s) """
    Cd = (0.49 + 0.065 * u10) * 1.e-3 # NOTE: there is an error in the paper. It has 0.0065 instead of 0.065
    Cd = np.where(u10<=11,0.0012,Cd)
    Cd = np.where(u10>=20,0.0018,Cd)
    return Cd

def calc_ustar(u10,rhosw=1026):
    """ Compute ustar from u10 (m/s) and density of sea water (rhosw in kg/m^3) """
    Cd = calc_cd10(u10)
    return np.sqrt(rhoair/rhosw * Cd) * u10

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 [3]:
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 [4]:
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

Compute friction velocity ($u_*$) and interpolate it and sea level pressure to CESM grid

$$ u_* = u_{10} \sqrt{\frac{\rho_a}{\rho_w} c_d} $$

$\rho_a$ = density of air (1.292 kg m$^{-3}$), $\rho_w$ = density of sea water, $c_d$ = drag coefficient, $u_{10}$ = wind speed at 10 m.

In [5]:
u10 = remap(u10_ncep,'u10','wind velocity at 10m','m/s')
slp = remap(slp_ncep.slp_,'slp','sea level pressure','Pa')
ustar = calc_ustar(u10,dsin.PD[:,0,...]) # compute ustar from u10
del ustar.coords['z_t']
ustar.attrs = {'long_name':'friction velocity', 'units':'m/s'}
ds = xr.Dataset({'ustar':ustar,'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 [6]:
fig, axs = plt.subplots(nrows=2,ncols=2,subplot_kw={'projection':proj},figsize=(9.5,6.3))
fig.subplots_adjust(wspace=0.05,hspace=0.25)
axs = axs.ravel()
varlist = ['TEMP','SALT','ustar','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),fontsize=10)
    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 [7]:
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 [8]:
fig, axs = plt.subplots(ncols=2,subplot_kw={'projection':proj},figsize=(9.5,4))
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',fontsize=10)
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}$)',fontsize=10)
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 [9]:
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=(9.5,4))
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$)',fontsize=10)
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.RdBu_r)
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$',fontsize=10)
cb = fig.colorbar(pm,ax=axs[1],orientation='horizontal',pad=0.03)
cb.set_label(r'%')

$^3\!$He gas fluxes

The total or net gas flux across the air-sea interface (mol m$^{-2}$ s$^{-1}$) is given by:

$$ F_{tot} = F_D + F_P + F_C $$

$F_D$ = diffusive gas exchange flux, $F_P$ = flux due to partially dissolved bubbles and $F_C$ = flux due to completely dissolved bubbles.

Diffusive gas exchange

The flux due to diffusion across the air-sea interface (mol m$^{-2}$ s$^{-1}$) is given by:

$$ F_D = K_s C_{eq} \left(\frac{P_{slp}}{P_{atm}} - \frac{C}{C_{eq}} \right) $$

$P_{slp}$ = local sea level pressure, $P_{atm}$ = standard atmospheric pressure, $C$ = $^3\!$He concentration at ocean surface.

$K_s$ is the gas transfer velocity across the ocean surface (m s$^{-1}$):

$$ K_s = \frac{u_*}{r_{wt} + \alpha\,r_{at}} $$

$r_{wt}$ and $r_{at}$ are the water-side and air-side resistance to transfer, respectively:

$$ r_{wt} = \sqrt{\frac{\rho_w}{\rho_a}} \left(h_w \sqrt{\frac{S_C}{\alpha_D}} + \frac{1}{\kappa}\,\ln{\frac{z_w}{\delta_w}} \right)$$$$ r_{at} = h_a \sqrt{S_{C_a}} + \frac{1}{\sqrt{c_d}} - 5 + \frac{\ln{S_{C_a}}}{2\kappa} $$

$\alpha$ is dimensioless solubility:

$$ \alpha = \frac{C_{eq}\,R\,T}{P_{atm}} $$

$S_C$ = Schmidt number of gas, $S_{C_a}$ = Schmidt number of air (0.9), $z_w$ = reference water depth (0.5), $\delta_w$ = molecular sublayer thickness (0.01), $\kappa$ = 0.4, $h_a$ = 0.13, $A$ = 1.3, $\phi$ = 1, $h_w = \frac{13.3}{A\,\phi}$, $R$ = ideal gas constant (8.314425 $\frac{\textrm{Pa m}^3}{\textrm{K mol}}$), $T$ = potential temperature (Kelvin). The $^4\!$He Schmidt number ($S_C$) is scaled by the $^3\!$He/$^4\!$He diffusivity ratio ($\alpha_D$) from Bourg & Sposito (2008).

In [10]:
sst    = ds.TEMP.isel(z_t=0)
rhow   = ds.PD.isel(z_t=0)
ustar  = ds.ustar
Pslp   = ds.slp
C      = ds.HELIUM3.isel(z_t=0)
Ceq    = ds.he3_ceq 
cd     = calc_cd10(ustar)

Sca    = 0.9
zw     = 0.5
deltaw = 0.01
kappa  = 0.4
ha     = 0.13
A      = 1.3
phi    = 1
hw     = 13.3 / (A * phi)

Sc = comp_he_schmidt(sst) / alpha_diff   # scale 4He Schmidt number by 3He/4He diffusivity ratio
rw = np.sqrt(rhow/rhoair) * (hw * np.sqrt(Sc) + np.log(zw/deltaw)/kappa)
ra = ha * np.sqrt(Sca) + 1/np.sqrt(cd) - 5 + np.log(Sca)/(2*kappa)
alpha = Ceq / Patm * R * (sst + T0K)
Ks = ustar / (rw + alpha * ra)
Fd = Ks * Ceq * (Pslp/Patm - C/Ceq)
Fd.attrs = {'long_name':'diffusive gas exchange','units':'mol/m^2/s'}
ds = ds.merge({'Fd':Fd})

Plot map

In [11]:
fig, (ax1, ax2) = plt.subplots(ncols=2,subplot_kw={'projection':proj},figsize=(9.5,4))
fig.subplots_adjust(left=0.05,wspace=0.025)
for ax in [ax1, ax2]:
    ax.set_global()
    ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')

Ks_cms = da2ma(Ks[0,...]) * 360000
pm = ax1.pcolormesh(lonc,latc,Ks_cms,transform=ccrs.PlateCarree())#,cmap=plt.cm.RdBu_r)
ax1.set_title(r'$K_s$',fontsize=10)
cb = fig.colorbar(pm,ax=ax1,orientation='horizontal',pad=0.03)
cb.set_label(r'cm h$^{-1}$')

Fd_pmol = da2ma(Fd.isel(time=0)) * spd * 1.e+12 # mol/m^2/s -> pmol/m^2/d
pm = ax2.pcolormesh(lonc,latc,Fd_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.RdBu_r)
center_cmap(pm)
ax2.set_title(r'Diffusive gas exchange ($F_D$)',fontsize=10)
cl = ax2.contour(loncp,ds.TLAT,Fd_pmol,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
cb = fig.colorbar(pm,ax=ax2,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))

min = -0.0128, max = 0.0138

Flux due to partially dissolved bubbles

The flux due to partially dissolved bubbles is given by:

$$ F_P = K_b C_{eq} \left( \left(1 + \Delta_P\right) \frac{P_{slp}}{P_{atm}} - \frac{C}{C_{eq}}\right) $$

$K_b$ is the gas transfer velocity (m s$^{-1}$) through large bubbles:

$$ K_b = 1.98 \times 10^6\, u_*^{2.76} \left(\frac{\frac{1}{\alpha_D}S_C}{660}\right)^{-\frac{2}{3}} \left[\frac{\textrm{cm}}{\textrm{hr}}\right] \times \frac{1}{360000} \left[\frac{\textrm{m}}{\textrm{cm}} \frac{\textrm{hr}}{s}\right] $$

The $^4\!$He Schmidt number ($S_c$) is scaled by the $^3\!$He/$^4\!$He diffusivity ratio ($\alpha_D$) from Bourg & Sposito (2008).

$\Delta_P$ is the supersaturation (as a fraction) due to partially dissolved bubbles:

$$ \Delta_P = \frac{152.44}{100}\,u_*^{1.06} $$

In [12]:
Kb = 1.98e+6 * np.power(ustar,2.76) * np.power(Sc/660,-2/3) * 1.e-2/3600 # cm/hr -> m/s
deltaP = 152.55 * np.power(ustar,1.06) / 100                            # % -> fraction
Fp = Kb * Ceq * ((1 + deltaP) * Pslp/Patm - C/Ceq)
Fp.attrs = {'long_name':'flux due to partially dissolved bubbles','units':'mol/m^2/s'}
ds = ds.merge({'Fp':Fp})

Plot map

In [13]:
fig, (ax1, ax2) = plt.subplots(ncols=2,subplot_kw={'projection':proj},figsize=(9.5,4))
fig.subplots_adjust(left=0.05,wspace=0.025)
for ax in [ax1, ax2]:
    ax.set_global()
    ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')

Kb_cms = da2ma(Kb[0,...]) * 360000
pm = ax1.pcolormesh(lonc,latc,Kb_cms,transform=ccrs.PlateCarree())
ax1.set_title(r'$K_b$',fontsize=10)
cb = fig.colorbar(pm,ax=ax1,orientation='horizontal',pad=0.03)
cb.set_label(r'cm h$^{-1}$')

Fp_pmol = da2ma(Fp.isel(time=0)) * spd * 1.e+12 # mol/m^2/s -> pmol/m^2/d
pm = ax2.pcolormesh(lonc,latc,Fp_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.RdBu_r)
center_cmap(pm)
ax2.set_title(r'flux due to partially dissolved bubbles ($F_P$)',fontsize=10)
cl = ax2.contour(loncp,ds.TLAT,Fp_pmol,[0],colors='k',transform=ccrs.PlateCarree(),linewidths=0.75)
cb = fig.colorbar(pm,ax=ax2,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))

min = -0.126, max = 1.08
In [14]:
fig, ax = plt.subplots(subplot_kw={'projection':proj})
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
data = da2ma(deltaP.isel(time=0)) * 100 # fraction -> %
ax.set_title(r'supersaturation due to partially dissolved bubbles ($\Delta_P$)',fontsize=10)
pm = ax.pcolormesh(lonc,latc,data,transform=ccrs.PlateCarree())
cb = fig.colorbar(pm,ax=ax,orientation='horizontal',pad=0.03)
cb.set_label(r'%')

Flux due to completely dissolved bubbles

The flux due to completely dissolved bubbles (mol m$^{-2}$ s$^{-1}$) is given by:

$$ F_C = 5.56\,u_*^{3.86}\,X_{He}\,I_r $$
In [15]:
Fc = 5.56 * np.power(ustar,3.86) * Xhe * Ir
Fc.attrs = {'long_name':'flux due to completely dissolved bubbles','units':'mol/m^2/s'}
ds = ds.merge({'Fc':Fc})

Plot map

In [16]:
fig, ax = plt.subplots(subplot_kw={'projection':proj})
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$)',fontsize=10)
pm = ax.pcolormesh(lonc,latc,Fc_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.RdBu_r)
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))

min = 4.22e-09, max = 0.484

Supersaturation at which the fluxes due to bubbles are balanced by the diffusive gas exchange

$$ \Delta_{eq} = \frac{K_b C_{eq} \Delta_P \frac{P_{slp}}{P_{atm}} + F_C}{(K_b + K_s) C_{eq} \frac{P_{slp}}{P_{atm}}} $$
In [17]:
delta_eq = (Kb * Ceq * deltaP * Pslp/Patm + Fc) / ((Kb + Ks) * Ceq * Pslp/Patm)

Plot map

In [18]:
fig, ax = plt.subplots(subplot_kw={'projection':proj})
ax.set_global()
ax.add_feature(cfeature.LAND,facecolor='#b0b0b0')
data = da2ma(delta_eq.isel(time=0)) * 100 # fraction -> percentage
ax.set_title(r'equilibrium supersaturation ($\Delta_{eq})$',fontsize=10)
pm = ax.pcolormesh(lonc,latc,data,transform=ccrs.PlateCarree())
cb = fig.colorbar(pm,ax=ax,orientation='horizontal',pad=0.03)
cb.set_label(r'%')

Total gas flux

In [19]:
Ft = Fd + Fp + Fc
fig, ax = plt.subplots(subplot_kw={'projection':proj})
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}$)',fontsize=10)
pm = ax.pcolormesh(lonc,latc,Ft_pmol,transform=ccrs.PlateCarree(),cmap=plt.cm.RdBu_r)
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))

min = -0.0669, max = 1.52
In [20]:
# ds.coords['time'].attrs = {} # something fishy going on here when I try to write the file
# ds.to_netcdf('testfoo.nc',encoding={'time':{'calendar':'noleap'}})

$K_s$ & $K_b$ as function of wind speed ($u_{10}$)

In [21]:
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(9.5,5))
sp = ax1.scatter(u10.values.ravel(),Ks.values.ravel()*360000,3,c=sst,marker='.') # m/s -> cm/hr
cb = fig.colorbar(sp,ax=ax1,orientation='horizontal')
_ = cb.set_label(r'SST ($^{\circ}$C)')
_ = ax1.set(xlabel=r'$u_{10}$ (m/s)',ylabel=r'$K_s$ (cm/hr)')
sp = ax2.scatter(u10.values.ravel(),Kb.values.ravel()*360000,3,c=sst,marker='.') # m/s -> cm/hr
cb = fig.colorbar(sp,ax=ax2,orientation='horizontal')
_ = cb.set_label(r'SST ($^{\circ}$C)')
_ = ax2.set(xlabel=r'$u_{10}$ (m/s)',ylabel=r'$K_b$ (cm/hr)')

$K_s$ & $K_b$ as function of wind speed ($u_{10}$) for $S_C = 660$

In [22]:
Sc = 660
rw = np.sqrt(rhow/rhoair) * (hw * np.sqrt(Sc) + np.log(zw/deltaw)/kappa)
ra = ha * np.sqrt(Sca) + 1/np.sqrt(cd) - 5 + np.log(Sca)/(2*kappa)
alpha = Ceq / Patm * R * (sst + T0K)
Ks = ustar / (rw + alpha * ra)
Kb = 1.98e+6 * np.power(ustar,2.76) * np.power(Sc/660,-2/3) * 1.e-2/3600 # cm/hr -> m/s
In [23]:
fig, axs = plt.subplots(nrows=2,ncols=2,figsize=(9.5,7.125))
fig.delaxes(axs[1,1])
axs= axs.ravel()
sp = axs[0].scatter(u10.values.ravel(),Ks.values.ravel()*360000,3,c='C0',marker='.') # m/s -> cm/hr
_ = axs[0].set(ylabel=r'$K_s$ (cm/hr)')
sp = axs[1].scatter(u10.values.ravel(),Kb.values.ravel()*360000,3,c='C0',marker='.') # m/s -> cm/hr
_ = axs[1].set(xlabel=r'$u_{10}$ (m/s)',ylabel=r'$K_b$ (cm/hr)')
Ktot = (Ks.values.ravel() + Kb.values.ravel()) * 360000
sp = axs[2].scatter(u10.values.ravel(),Ktot,3,c='C0',marker='.') # m/s -> cm/hr
_ = axs[2].set(xlabel=r'$u_{10}$ (m/s)',ylabel=r'$K_s + K_b$ (cm/hr)')