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

Soil properties

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
  • Soil texture ($tex$):
    • 0 if clay
    • 9 if sand
    • 0.07 if area type is land_ice
  • 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.

Albedo Calculation

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.

Bare soil

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$$

Land Ice

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}$$

Surface drag parameters

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

Roughness length

$$ z_0 = \left \{ \begin{array}{1} 0.01m,\ (soil) \\ 0.001m,\ (snow\ or\ ice)\\ 0.0002m,\ (sea) \end{array} \right .$$

Neutral drag coefficient

$$ C_{DN} = \left ( \frac{\kappa}{\log{z_m} - \log{z_0}} \right )^2 ,$$

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

Normalised mixing length

$$\zeta = \exp{\left (\frac{-\kappa}{\sqrt{C_{DN}}}\right )}$$

Richardson number

$$Ri = -\frac{gz_m}{T_s U_m^2}\left [ T_u - T_s\right ] $$

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.

Momentum drag coefficient

$$ C_{Dm} = \left \{ \begin{array}{1} C_{DN}\left [ 1 - \frac{8Ri}{1 + 56.768 C_{DN}\sqrt{\frac{-Ri}{\zeta}}}\right ],\ Ri<0 \\ C_{DN} \left [\frac{(1 - 4\epsilon Ri)^2}{1 + 8(1-\epsilon)Ri} \right ],\ Ri \ge 0 \end{array} \right .$$

Heat and other scalars drag coefficient

$$ C_{Dh} = \left \{ \begin{array}{1} C_{DN}\left [ 1 - \frac{12Ri}{1 + 41.801 C_{DN}\sqrt{\frac{-Ri}{\zeta}}}\right ],\ Ri<0 \\ C_{DN} \left [\frac{1 - 4\epsilon Ri}{1 + (6-4\epsilon)Ri} \right ]^2,\ Ri \ge 0 \end{array} \right .$$

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

Surface fluxes

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

Sensible heat flux (Eq. 8.3)

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.

Latent heat flux (Eq. 8.4)

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}$$


$$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 )$$


  • $\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.

Momentum fluxes

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

Surface energy balance

The energy balance is straightforward once the latent and sensible heat fluxes are known. Adding the radiative fluxes completes the balance.

$$Net = Net_{SW} + Net_{LW} - LHF - SHF$$

Moisture and heat conduction through snow, ice and soil

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$$


$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.

Bottom Boundary conditions

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

Numerical scheme

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.