`SecondBEST`

is a simplified implementation of the Bare Essentials for Surface Transfer (BEST) land surface model. The main simplifications are that `SecondBEST`

does not incorporate the effects of vegetation, which reduces a lot of calculations within BEST. We also assume each grid cell is either fully bare soil or snow/ice, i.e, there are no fractional land types. The implementation here follows the description in Pitman et al.

This notebook is a ready reference to the equations used in `SecondBEST`

so that interested users don't have to read the above reference if they don't wish to. However, if you plan on making heavy use of `SecondBEST`

, you are encouraged to do so!

**Surface types**:

*s*= air above the surface, or lowest atmospheric model layer.*u*= upper soil layer. BEST uses three soild layers.*l*= lower soil layer.*b*= base soil layer.*n*= snow.*d*= intercepted water.

**Spectral regions, phases of water**

*L*= liquid.*F*= frozen.*SW*= shortwave.*LW*= longwave.

**Water content**

- $X_w$ = volumetric soil liquid water content
- $X_i$ = volumetric soil ice content

Certain soil properties are calculated initially based on the soil and area types.
Currently, `SecondBEST`

recognises two soil types: `clay`

and `sand`

.

- Soil colour ($clr$):
- 0.2 if
`clay`

- 1.0 if
`sand`

- 0.2 if
- Soil texture ($tex$):
- 0 if
`clay`

- 9 if
`sand`

- 0.07 if area type is
`land_ice`

- 0 if
- Soil porosity ($X_v$), Eq. 4.10:
- $0.6 - 0.03\times tex$

- Soil moisture volume fraction at field capacity ($X_{FC}$), Eq. 4.11:
- $(0.95 - 0.086\times tex)X_v$

- Soil moisture volume fraction at wilting point ($X_{wilt}$), Eq. 4.12:
- $X_v - 0.03$
- 0.01 when area type is
`land_ice`

.

The albedo calculation is described in section 5.0 of Pitman et al. `SecondBEST`

calculates albedo in the near-infraread and shortwave separately. This is also quite convenient since RRTMG, the radiative transfer code also uses two separate albedos.

The albedo for bare soil is calculated as (Eqns. 5.5, 5.6): $$\alpha_{SW,u} = 0.10 + 0.1clr + 0.06(1-W_{L,u})\\ \alpha_{LW,u} = 2\alpha_{SW,u}$$

where $clr$ is the soil colour defined above.

$W_{L,u}$ is the soil wetness factor, and is defined as the ratio of the liquid soil water content in the upper soil layer to its porosity.

$$W_{L,u} = X_{w,u}/X_v$$The albedo of Land Ice is calculated as (Eqns. 5.7, 5.8): $$\alpha_{SW,u} = 0.60 + 0.06(1-W_{L,u}) \\ \alpha_{LW,u} = \frac{1}{3}\alpha_{SW,u}$$

`SecondBEST`

calculates the drag parameters as in section 6 of Pitman et al.

where $z_m$ is the hydrostatic mid-level height of the lowest atmosphere level and $\kappa$ is the von Karman constant.

Where $U_m$ is the magnitude of the wind speed at the lowest model level and the other notation is as before. $U_m$ has a minimum value that can be set during initialisation of `SecondBEST`

.
A negative Richardson number indicates unstable conditions.

where $\epsilon$ is 0.01 for land and 1.0 for oceans.

The equations are from Section 8. in Pitman et al.

is given by the relatively well known bulk aerodynamic formula.

$$SHF = \rho_s C_{pd} U_m C_{Dh} (T_u - T_s)$$where $C_{pd}$ is the heat capacity at constant pressure of dry air, and $\rho_s$ is the density of air at the surface.

is also given by the bulk formula

$$LHF = L_v\rho_s U_m C_{Dh} \beta_u(q^*_u - q_s) $$where $q^*_u$ is the saturation specific humidity at the soil surface. A "wetness factor", $\beta_u$, is also required, defined as

$$\beta_u = \left \{ \begin{array}{1} 1,\ snow\\ \frac{W_{F,u}L_v}{L_i} + \frac{E_{usmax}}{E_{usP}},\ soil\end{array} \right .$$$L_v$ is the latent heat of vapourisation, $L_i$ is the latent heat of sublimation, $W_{F,u}$ is the frozen soil moisture in the upper soil layer. The first term containing these terms represents the rate at which frozen soil can sublimate.

$E_{usP}$ is the potential evaporation rate at the soil surface,

$$E_{usP} = \rho_s c_u (q^*_u - q_s)$$where $c_u$ is the soil conductance, defined as (Eq. 7.22)

$$c_u = C_{Dh}U_m.$$$E_{usmax}$ is the exfiltration rate given by

$$E_{usmax} = K_{HD}\Theta_u^{0.5B + 2} - K_{H0}\Theta_u^{2B + 3}$$Where

$$B = \left \{ \begin{array} 10,\ clay\\ 4,\ sand\end{array} \right .$$$$K_{H0} = \left \{ \begin{array} 0.001,\ clay\\ 0.1,\ sand \end{array} \right .$$and $$K_{HD} = \left ( \frac{-4K_{H0}B\Psi_0\rho_wX_v(1 - W_{F,u})}{\pi\Delta t}\right )$$

where

- $\rho_w$ is the density of the water
- $\Psi_0 = -0.2m$ for all soil types is the soil suction at saturation.
- $\Delta t$ is the model time step
- $\pi$ is the constant.

$\Theta_u$ is defined as

$$\Theta_u = \frac{W_{L,u} - 0.01}{1 - W_{F,u}}$$The wetness factors $W_{L,x}$ in the soil in the upper and lower layers is **NOT** allowed to fall below 0.01 due to physical and numerical issues.

are directly calculated from the momentum bulk coefficient $$U_s = -\rho_s C_{Dm}U_mu$$

The transfer of heat and moisture through the snow-ice-soil pack requires solving the diffusion equation along with local sources and sinks of heat and moisture which correspond to freezing/melting of water.

The conservation laws are:

$$c_vT_t = G_z + L_f\Gamma + c_{pw}RT_z - T_t(fict)$$$$\rho_w (X_w)_t = (R - (1-f_i)E)_z - \Gamma - \Pi$$$$\rho_i(X_i)_t = -(f_i E)_z + \Gamma$$where

$f_i = X_i/(X_w + X_i)$, $c_v$ is a volumetric heat capacity for soil.

$\Gamma = \frac{c_v}{L_f}\frac{(T_f - T)}{\Delta t}$ is the potential ice production rate, which is zero if there is no water for freezing or ice for melting.

$\Pi = A_vE_c/\Delta z$ is the water removed by transpiration, which is zero in this component (no vegetation effects)

And the flux densities are

$G = K_T T_z$

$R = \rho_w K_H(1 + \psi_z) = \rho_w K_H(z + \psi)_z$ for water

$E = \left \{ \begin{array}\\ 0,\ z < 0\\ E_{us},\ z=0 \end{array} \right .$ for water vapour

$T_t(fict)$ is a fictional source that is not necessary in our case since we are solving the full heat/mass transfer equations.

We use a Neumann condition for all variables $()_z=0$, representing a sink of energy and water to the deep soil.

We use an implicit scheme to step forward the diffusive (density) terms and an Euler scheme to step forward the other sources and sinks. The term $c_{pw}RT_z$ represents advection, and this will set the CFL condition $\Delta z/(c_{pw}R) > \Delta t/2$ where $\Delta t$ is the timestep to be used.