PlaSim/Plasim_RM_16/simba.tex

\chapter{Dynamic Vegetation}\label{dynveg}
A simple terrestrial dynamic global vegetation model (DGVM), Simulator
for Biospheric Aspects (SimBA), is used to obtain the following land
surface variables for non-glaciated grid cells: surface albedo $A$,
the roughness length $z_{0}$, a surface conductance factor for the
latent heat flux $C_{w}$, and a "bucket" depth for the soil, $W_{max}$.
These land surface variables depend on SimBA variables-- the latter of
which ultimately depend on the following three  (global) variables:
soil moisture content ($W_{soil}$), snow depth ($z_{snow}$), and
vegetative biomass ($BM$).  Of these three variables, vegetative
biomass has the greatest importance within SimBA.

\section{Equations for SimBA Variables}
Vegetative biomass ($BM$) depends on net primary productivity as given in SimBA's fundamental governing equation:

\begin{equation}
\label{eq:dBMdt} 
\frac{\partial BM}{\partial t} = NPP - \frac{BM}{\tau_{veg}}
\end{equation}
where $\tau_{veg}$ is the residence time of the vegetative carbon and equals 10 years, and $NPP$ (net primary productivity) is approximated as $0.5 * GPP$.  

The approximation $NPP = 0.5 * GPP$ is briefly justified in \cite{kleidon2006}, but  some recent studies find that $NPP/GPP$ can deviate considerably from 0.5 (\cite{delucia2007}; \cite{zhang2009}).  Nevertheless, the $NPP/GPP$ = constant parameterization is attractively simple, and it has been assumed by widely-used productivity models such as CASA and FOREST-BGC (\cite{delucia2007}).  A scheme for heterotrophic respiration is currently included in the model as a diagnostic only.  The gross primary production formulation is detailed in the next subsection.

\subsection{Gross Primary Production}

$GPP$ is calculated as the minimum of a water-limited rate and a light-limited rate.  (That is, $GPP$ = $min(GPP_{light},GPP_{water}$.)  This approach originates in a crop model (\cite{monteith1989}) which was later adapted for forest canopies (\cite{dewar1997}).  

\subsubsection{Light-limited Gross Primary Production}

The light-limited rate, $GPP_{light}$ follows a light-use efficiency approach (e.g. see \cite{yuan2007}) as follows:

\begin{equation}
\label{eq:gppl}
GPP_{light} = \epsilon_{luemax} * \beta(CO_{2}) * f(T_{sfc}) * fPAR * SW\downarrow
\end{equation}
where $\epsilon_{luemax}$ is a globally-constant maximum light use efficiency parameter = $3.4*10^{-10}$~kgC/J,
$\beta(CO_{2})$ represents a logarithm-based "beta" factor effect on productivity for when $CO_2$ concentration deviates from the reference value of 360ppmv (see below),
$f(T_{sfc})$ is a temperature limitation function (defined below) which lowers productivity for cold temperatures, 
$fPAR$ is the fraction of photosynthetically active radiation that is absorbed by green vegetation (see below), and $SW\downarrow$ is the downward flux of shortwave radiation at the surface (in $W/{m^2}$).

In equation~\eqref{eq:gppl}, the first term on the right hand side, $\epsilon_{luemax}$, is the light use efficiency with respect to the absorbed total shortwave broadband radiation.  The value of $3.4*10^{-10}$~kgC per J (of fPAR * $SW\downarrow$) is derived from the maximum light use efficiency value of the CASA model, 0.389 gC MJ$^{-1}$ of APAR [absorbed photosynthetically active radiation] (\cite{potter1993};\cite{field1995}) by using the commonly-used approximations GPP = 0.5*NPP and $SW\downarrow$ = 0.5 * PAR [photosynthetically active radiation] (at top of the canopy).  (Both approximations are made in the CASA model (\cite{potter1993};\cite{field1995})).  The equivalent $\epsilon_{luemax}$ value in SimBA would be $3.89*10^{-10}$~kgC/J, but this is lowered to $3.4*10^{-10}$~kgC/J to account for the lack of an optimum growing temperature in SimBA, since the lack of such an optimum causes light-limited productivity to be slightly overestimated for most regions.

The second term in equation~\eqref{eq:gppl}, $\beta(CO_{2})$, is taken from (\cite{harvey1989}), but incorporates carbon compensation point as follows:

\begin{equation}
\label{eq:betafactor}
\beta(CO_{2}) = 1 + \max{(0,BF * \ln{(\frac{CO_{2}-CO_{2,comp}}{CO_{2,ref}-CO_{2,comp}}}))}
\end{equation}
where BF = the carbon dioxide sensitivity or "beta" factor, $CO_{2,ref}$ = 360ppmv, and $CO_{2,comp}$ = the light compensation point (in ppmv) (set to zero by default).

The third term in equation~\eqref{eq:gppl}, $f(T_{sfc})$, is as follows:

\begin{equation}
f(T_{sfc}) = \left\{ \begin{array}{ll} 
  0 & \mbox{if $T_{sfc} \leq 0\,^{\circ}\mathrm{C}$} \\
  \frac{T_{sfc}}{T_{crit}} & \mbox{if $ 0<T_{sfc}<T_{crit}$} \\
  1 & \mbox{if $T_{sfc} \geq T_{crit}$}
\end{array} \right.
\end{equation}
where $T_{sfc}$ is the surface temperature in $^{\circ}\mathrm{C}$, and $T_{crit}$ is the "critical" temperature (set to 5$\,^{\circ}\mathrm{C}$) by default) at which temperature is no longer limiting to productivity.

The fourth term in equation ~\eqref{eq:gppl}, fPAR, is also referred to as "vegetation cover" ($f_{veg}$) in the model.  ``fPAR'' refers to fraction of photosynthetically active radiation (PAR) that is absorbed by photosynthesizing parts (i.e. green leaves) of plants.  (PAR is the fraction of incoming solar radiation that is in wavelengths usable for photosynthesis.)  In a Beer's Law approach, fPAR can be approximated as a function of the leaf area index (LAI):

\begin{equation}
\label{eq:fPAR}
fPAR = 1 - e^{-k_{veg} * LAI}
\end{equation}
where $k_{veg}$ is a light extinction coefficient (set to 0.5 by default).  This same approach is followed in the forest canopy and crop models on which SimBA is based (\cite{monteith1989}; \cite{dewar1997}), and it is also used in the formulation of snow-free surface albedo in the ECHAM5 GCM (\cite{rechid2009}) as well as in SimBA (see section below on surface albedo).

\subsubsection{Water-limited Gross Primary Production}

The water-limited rate, $GPP_{water}$, whose equation we derive here, is based on the equivalent formulation in the forest canopy model of \cite{dewar1997} and follows a ``big leaf"-diffusivity approach.  Diffusions of CO$_{2}$ and H$_{2}$O between leaf and atmosphere are proportional to the concentration gradient between leaf intercellular air and atmosphere times the respective diffusivities of these gases.  In our ``big leaf'' model, the canopy is treated as if it were a large single leaf that is well coupled to the atmosphere.  Such strong leaf-atmosphere coupling permits us to neglect the leaf boundary layer conductance as compared to the stomatal conductance.  Following the appendix section of \cite{dewar1997}, we can now write an equation for water-use efficiency at the leaf scale, $q_{leaf}$:

\begin{equation}
\label{eq:qleaf}
q_{leaf} = \frac{c_{a}}{1.6VPD}(1-\frac{c_{i}}{c_{a}})
\end{equation}
where ${c_{i}}$ is the intercellular CO$_2$ concentration, ${c_{a}}$ is the atmospheric CO$_2$ concentration, VPD = the vapor pressure deficit between the (presumed) saturated leaf surface and atmosphere, and the ``1.6'' term represents the difference in diffusivity between CO$_{2}$ and H$_{2}$O due to their differing molecular weights.  Here,  the units are assigned as follows: $q_{leaf}$ is in $\frac{mol C}{mol H_{2}O}$; ${c_{i}}$ and ${c_{a}}$ are in Pa CO$_2$; VPD is in Pa H$_2$O; and the ``1.6'' term has units of $\frac{\frac{mol H_{2}O m^{-2} s^{-1}}{Pa H_{2}O/Pa air}}{\frac{mol CO_{2} m^{-2} s^{-1}}{Pa CO_{2}/Pa air}} = \frac{Pa CO_2}{Pa H_{2}O} \frac{mol H_{2}O}{mol CO_2} = \frac{Pa CO_2}{Pa H_{2}O} \frac{mol H_{2}O}{mol C}$ (where ``Pa air" is pascals of total atmospheric pressure, and where the last equality holds because there is one mole of carbon in each mole of CO$_2$).

Equation~\eqref{eq:qleaf} can be written as mass- rather than mole-based water-use efficiency by multiplying the right hand side by a conversion factor of $\frac{12 kg C / 1000 mol C}{18 kg H_2O / 1000 mol H_2O}$, i.e. by $\frac{2}{3} \frac{kg C mol H_2O}{kg H_2O mol C}$.  We call the new left hand side ``$q_{mass}$'' (with units of $\frac{kg C}{kg H_2O}$):

\begin{equation}
\label{eq:qmass}
q_{mass} = \frac{c_{a}}{1.6VPD}(1-\frac{c_{i}}{c_{a}})
\end{equation}

Next, we make a simplifying assumption that $\frac{c_i}{c_a}$ has a constant value of 0.7 to get

\begin{equation}
\label{eq:qmass2}
q_{mass} = \frac{2}{3}\frac{0.3 c_{a}}{1.6VPD}
\end{equation}

It should be noted that the simplifying  assumption that
$\frac{c_i}{c_a}$ is constant is not entirely valid.
Although early studies showed that the ratio is largely conserved
as ambient CO$_2$ is varied (e.g. \cite{wong1979}; \cite{polley1993}),
many recent studies have found significant variance of the ratio
(e.g. from 0.54 to 0.95 (\cite{brooks1997}), such that decreasing
$\frac{c_i}{c_a}$ is associated with increasing light (\cite{brooks1997}),
increasing leaf-atmosphere vapor pressure deficit (\cite{morison1983}),
greater soil moisture stress (\cite{turnbull2002}), and higher canopy
position (\cite{brooks1997}).  In addition, C$_4$ plants tend to have
lower values than C$_3$ plants (e.g. see \cite{bunce2005}), which should
thus cause SimBA to underestimate productivity in climatic zones which
favor C$_4$ plants (warm and dry areas).  On the other hand, at least
one other more sophisticated productivity model assumes a constant
$\frac{c_i}{c_a}$ (\cite{knorr2000}). 

We can now use the relation, $NPP_{water} = 0.5 \ GPP_{water}$,
and $q_{mass}$, water-use efficiency, from equation~\eqref{eq:qmass2} to write a large scale
expression for water-limited gross primary production ($GPP_{water}$) as follows:

\begin{equation}
\label{eq:gppwbasic}
GPP_{water} = 2 \ q_{mass} \ T
\end{equation}
where T = transpiration (in $kg H_2O \ m^{-2} s^{-1}$), and $GPP_{water}$ has units of $kg C m^{-2} s^{-1}$.

Next, we assume, as is also done in the original MOSES land surface model for the no-wet canopy case (\cite{cox1999}), that transpiration's contribution to total evapotranspiration (ET) equals the vegetative cover fraction, $f_{veg}$, i.e. that

\begin{equation}
\label{eq:transpiration}
T = ET * f_{veg}
\end{equation}
(Recall that we define $f_{veg} = fPAR = 1 - e^{- k{veg} * LAI}$ in SimBA.)

Next, we substitute in equations ~\eqref{eq:qmass2} and ~\eqref{eq:transpiration} into equation~\eqref{eq:gppwbasic}; use the conversion $c_a = P * co2 * 10^{-6}$, where P and $c_a$ are atmospheric pressure and atmospheric carbon dioxide partial pressure, respectively (in Pa), and where co2 is the atmospheric carbon dioxide concentration (in ppmv); and use a reference density for H$_2$O of 1000 kg m$^{-3}$ to get our final form equation for water-limited gross primary production:

\begin{equation}
\label{eq:gppw}
GPP_{water} = \frac{co2conv * P * f_{veg} * ET * (0.3 * co2)}{VPD}
\end{equation}
where co2conv = 8.3 * 10$^{-4}$, ET is evapotranspiration (in $\frac{m^3}{m^2 s}$), and GPP$_{water}$ is again in $kg C m^{-2} s^{-1}$.

The above GPP$_{light}$ and GPP$_{water}$ formulas rely on vegetative cover fraction (f$_{veg}$) (which equals fPAR), and f$_{veg}$ is a function of LAI.  The next subsection describes how these are derived in SimBA.

\subsection{Vegetative Cover}

% PP: this subsection relies heavily on the old documentation

The vegetative cover $f_{veg}$ of the land surface, that is the fraction which
is covered by green biomass (leaves), is computed as the minimum of an 
water-limited value
$f_{veg,w}$ and a structurally limited value $f_{veg,s}$:

\begin{equation}
f_{veg} = \min{(f_{veg,w}, f_{veg,s})}
\end{equation}

The water-limited vegetation cover $f_{veg,w}$ is computed from the soil moisture $W_{soil}$ as:

\begin{equation}
f_{veg,w} = f_{W_{soil}}
\end{equation}
with the function $f_{W_{soil}}$ given by:

\begin{equation}
\label{eq:waterstressfactor}
f_{W_{soil}} = \min{(1, \max{(0, \frac{W_{soil} / \! W_{max}}{W_{crit}})})}
\end{equation}
where $W_{crit}$ = 0.25 and $W_{max}$ is the biomass-dependent soil "bucket" depth (see below for derivation).  This water limitation function, 
$f_{W_{soil}}$, is motivated by the fact that water stress for plants sets in at a critical 
value $W_{crit}$.  For simplicity, a fractional water content is used rather than 
a specific matric potential which would reflect the permanent wilting point.  Other land surface models use a similar approach (e.g. \cite{albertson2001})

The structurally limited vegetation cover $f_{veg,s}$ is obtained from a structurally limited maximum 
leaf area index $LAI_{m}$ as follows:

\begin{equation}
f_{veg,s} = 1 - \exp{(-k_{veg} \, LAI_{m})}
\end{equation}
which is sustained with the present amount of biomass:

\begin{equation}
LAI_{m} = 0.1 + \frac{2}{\pi} * LAI_{max} * atan(c_{{veg}_l} * BM)
\end{equation}
where LAI$_{max}$ represents the theoretical LAI that is approached as biomass becomes infinitely large, and with the values for LAI$_{max}$ = 9 and c$_{{veg}_l}$ = 0.25 giving the best fit to the old scheme used in earlier versions of SimBA for mid to high biomass (BM) values, while maintaing realism for low BM values.  (See notes at end of section on forest cover for more info on how LAI parameterization was developed.)

$LAI = -\frac{\ln{(1 - f_{veg})}}{k_{veg}}$ inputs back the effect of water limitation to LAI.  As a final note, in the scheme presented here, LAI follows drought-deciduous phenology, but it does not follow winter-deciduous phenology.

\subsection{Forest Cover}

Forest cover is an influential variable in the model.  In SimBA, forest cover refers to non-prostrate woody vegetation which sticks out above the snow pack.  Only forest cover contributes to surface roughness, and only forest cover lowers the snow-covered surface albedo relative to that of snow-covered bare soil.  Non-forested but vegetated land acts the same as bare soil with respect to surface roughness and snow-covered albedo.

Forest cover (F) is parameterized as follows:

\begin{equation}
F = 0.1 + \max(0,\frac{atan(BM - c_{{veg}_a}) - atan(c_{{veg}_d})}{0.5\pi - atan(c_{{veg}_d})})
\end{equation}
where 
$c_{{veg}_a}$ + $c_{{veg}_d}$ = $f_{crit}$, where $f_{crit}$ is defined as the biomass threshold at which forest cover begins to rise above zero.

$c_{{veg}_a}$ and $c_{{veg}_d}$ are currently set at 2.9 and -1.9, respectively, and thus $f_{crit}$ = 1.0 kg m$^{-2}$.  $c_{{veg}_a}$ is chosen to be near 3.0, which was the analogous value for the previous version of SimBA, thus giving an excellent fit between new and old forest cover parameterizations for biomass values above 3 kg m$^{-2}$.

\subsubsection{How $f_{crit}$ Was Chosen} 

Based on NPP and ecosystem type model data 
(\cite{mcguire1992}; \cite{cramer1999}), it was estimated that woody shrub cover begins to occur at circa 100 g per m$^2$ per year of NPP in the real world, which equals 1 kg of biomass per m$^2$ in Simba under steady state conditions.
This is a rough estimate, because forest cover is not solely a function of NPP.

\subsubsection{More on the Basis for the Biomass-LAI and Biomass-forest Cover Relationships}

The relationships were based more on NPP than on biomass.  NPP data from \cite{mcguire1992} and \cite{cramer1999} was converted into biomass data by using the steady state approximation $NPP \approx \frac{BM}{\tau_{veg}}$, which is obtained when the LHS of equation~\eqref{eq:dBMdt} is approximately 0 over long time scales.
Forest cover and LAI were related by using the land surface dataset for the ECHAM GCM (\cite{hagemann2002}).  Note that  \cite{hagemann2002} apparently does not consider woody shrub cover to be partial forest cover, as does SimBA, and this was taken into account in creating the two formulations.
Finally, LAI data for some arctic and mountain ecosystems (\cite{bliss1981}, p.195 and p.219) was also used to calibrate this part of the model.

\section{Derivation of Land Surface Parameters}

\subsection{Soil Water Holding Capacity} We informally refer to it here as soil ``bucket'' depth.  The general idea behind this formulation is that the bucket depth increases as the root biomass increases.  This dependence on root biomass has been incorporated into other simple land surface schemes (e.g. ENTS for the GENIE climate system model (\cite{williamson2006})).

The non-linear relationship between biomass and bucket depth (W$_{max}$) is due to the non-linear dependence of an intermediate variable, Vsoil, on biomass.  The Vsoil-biomass relationship was originally designed to be perfectly analogous  to the forest cover-biomass relationship, except for the former's using root biomass and the latter's using "shoot" biomass.  This updated version of SimBA maintains essentially the same Vsoil-biomass relationship, except that now it depends on total biomass.  (Note: the code can be easily modified to account for root-shoot partitioning different than the implicit 0.5 value.)  Here is the Vsoil-biomass equation:

\begin{equation}
Vsoil = \min(1.0,\max(0,\frac{atan(BM - c_{{veg}_f}) - atan(-c_{{veg}_f})}{0.5\pi - atan(-c_{{veg}_f})}))
\end{equation}
where $c_{{veg_f}}$ is set to 3.0 by default.

Next, soil bucket depth, W$_{max}$ depends linearly on Vsoil as follows:

 \begin{equation}
W_{max} = W_{max_{max}} * Vsoil + W_{max_{min}} * (1 - Vsoil)
\end{equation}
where W$_{max_{max}}$ = 0.5 is the theoretical soil bucket depth as biomass becomes infinitely large and W$_{max_{min}}$ = 0.05 is the soil bucket depth when biomass = 0 kg m$^{-2}$ (i.e. for bare soil).  Each bucket depth has a unit of meters in the model, and the given values for W$_{max_{max}}$ and W$_{max_{min}}$ are taken from \cite{kleidon2006}.

\subsection{Surface Albedo}

Surface albedo is first calculated for snow-free conditions, and then it is modified if there is snow.  Solar zenith angle dependence, diffuse-direct radiation distinction, and the dependence of bare soil albedo on soil moisture content are all neglected.

\subsubsection{Snow-free Surface Albedo}

This formulation is identical to that used in ECHAM5 (\cite{rechid2009}) and is as follows:

\begin{equation}
A_{snow-free} = A_{fully-leaved} * f_{veg} + A_{bare} * (1-f_{veg})
\end{equation}
where ``A'' denotes albedo; ``fully-leaved'' and ``bare'' denote conditions of infinite LAI and zero LAI, respectively; and f$_{veg}$ is vegetation cover (as before) = $1 - e^{-k_{veg}*LAI}$, where $k_{veg}$ is -0.5 as before.  A$_{fully-leaved}$ and A$_{bare}$ are currently set to 0.12 and 0.30, respectively.

It is important to note that snow-free albedo depends only on leaf area index (LAI) and not on forest cover.  Stems and branches, etc., are tacitly assumed to have the same albedo as bare soil, 0.3.

\subsubsection{Surface Albedo when Snow is Present} 

The grid cell is divided up into a forest-covered part and non-forest covered part.  Non-forest cover is a mixture of prostrate vegetation (e.g. grass, non-shrubby tundra) and bare soil.  The albedo of snow-free non-forest cover is mixed in with that of deep snow concurrent with the snow depth as follows:

\begin{equation}
A_{NF} = A_{NF_{snow-free}} + (A_{snow} - A_{NF_{snow-free}}) \frac{snowdepth}{snowdepth + 0.01}
\end{equation}
where ``NF'' denotes non-forest cover, A$_{snow}$ is the albedo of deep snow, and snowdepth is in meters.  Note:  the albedo of the snow-free non-forest cover is tacitly taken to equal the albedo of the entire grid cell under snow-free conditions. 

Deep snow albedo, A$_{snow}$, lowers with increasing surface temperature (T$_{sfc }$) as follows:

\begin{equation}
A_{snow} = \left\{ \begin{array}{ll} 
  A_{snow_{max}} & \mbox{if $T_{sfc} \leq -10\,^{\circ}\mathrm{C}$} \\
  A_{snow_{min}} + (A_{snow_{max}} - A_{snow_{min}}) (\frac{T_{sfc }}{-10\,^{\circ}\mathrm{C}}) & \mbox{if $ -10\,^{\circ}\mathrm{C} < T_{sfc }< 0\,^{\circ}\mathrm{C}$} \\
  A_{snow_{min}} & \mbox{if $T_{sfc} = 0\,^{\circ}\mathrm{C}$}
\end{array} \right.
\end{equation}
where A$_{snow_{max}}$ = 0.8 and A$_{snow_{min}}$ = 0.4.

Forest cover (again denoted by ``F'') is modeled to protrude from the snow pack and mask the snow beneath it.  For simplicity, the forest-covered portion of the grid cell is assigned the same albedo, A$_{F_{snow}}$, regardless of surface temperature or the amount of snow accumulation.  A$_{F_{snow}}$ is assigned a default value of 0.20 in the model.  Earlier versions of SimBA had A$_{F_{snow}}$ = 0.35.  The lower value of 0.20 has been adopted for a number of reasons:  it is used in the ECHAM GCM for fully snow-covered evergreen forests (\cite{roesch2001}; \cite{roesch2006}); and it is closer to satellite (\cite{gao2005}) and field (\cite{betts1997}) measurement values, particularly when the trees are not cold winter-deciduous (as is the case in SimBA).

Finally, the albedo ``A'' for the entire grid cell, is taken to be the linear combination of the respective albedos for forest-covered (F) and non-forest-covered (NF) fractions:

\begin{equation}
A = A_{F_{snow}} * F +  A_{NF} * (1 - F)
\end{equation}

\subsection{Surface Conductance}

Also denoted as ``rhs'' or as ``surface wetness'' within the model, we denote it here as ``C$_w$''.

\begin{equation}
C_w = pgs * f_{W_{soil}} 
\end{equation}
where f$_{W_{soil}}$ is taken from equation ~\eqref{eq:waterstressfactor} and represents a water stress factor due to reduced soil moisture content, and pgs represents the surface conductance achieved under non-water-stressed conditions.  The default value for pgs in the model is 1.0.  Finally, as a correction from previous versions of SimBA, C$_w$ is set to 1.0 regardless of f$_{W_{soil}}$ when snowcover $> 0$, to take into account the presence of sublimatable snow at the surface.

\subsection{Surface Roughness}

Surface roughness is taken as a non-linear combination of roughness due to orography and roughness due to vegetation.  As mentioned earlier, surface roughness due to vegetation is a function of forest cover only.  Hence, no increase in surface roughness occurs as biomass goes from 0 kg m$^{-2}$ to 1 kg m$^{-2}$ (the value at which forest cover commences).  We denote surface roughness due to vegetation as z$_ {0,veg}$ and formulate it as follows: 

\begin{equation}
z_{0,veg} = F * (z_{0,F}) + (1 - F) \, (z_{0,NF})
\end{equation}
where ``F'' and ``NF'' denote ``forest cover'' and ``non-forest cover'', respectively; $z_{0,NF}$=0.05~m, the vegetative surface roughness in the absence of forest cover; and $z_{0,F}$=2~m, the vegetative surface roughness when fully-forested (i.e. when F = 1).

Finally, surface roughness of a grid cell, z$_0$, is formulated as follows:

\begin{equation}
z_0 = \sqrt{{z_{0,veg}}^2 + {z_{0,oro}}^2}
\end{equation}
where z$_{0,oro}$ is the surface roughness due solely to orography.