Land Surface Model¶
WSIM includes a Land Surface Model (LSM), based on a synthesis of two frequently cited soil moisture/runoff models: WBM [2] and “Leaky Bucket” [12][25]. The LSM is driven by three variables:
Average monthly temperature (\(T\))
Total monthly precipitation (\(Pr\))
The fraction of days in the month during which precipitation falls (\(p_{wet}\))
It produces estimates of the following quantities:
Parameter 
Symbol 
Variable Name 

Evapotranspiration 
\(E\) 

Evapotranspiration minus potential evapotranspiration 
\(EE_0\) 

Net precipitation 
\(P\) 

Potential evapotranspiration 
\(E_0\) 

Potential minus actual evapotranspiration 
\(E_0E\) 

Runoff 
\(R\) 

Soil Moisture 
\(W_s\) 

Snow Accumulation 
\(S_a\) 

Snow Melt 
\(S_m\) 

Total Blue Water (runoff plus upstream runoff) 

Note
Runoff quanties (RO_mm
and Bt_RO
) are also output in variant
forms (Runoff_mm
and Bt_Runoff
) that represent runoff computed
for a model iteration in isolation, ignoring all effects of runoff
detention between timesteps.
Core functionality of the LSM is implemented in the wsim.lsm
R package.
Data dependencies¶
In addition to the forcing data (\(T\), \(Pr\), \(p_{wet}\)), the LSM relies on the following static data:
elevation, used to estimate the rate of snow melting
soil total available water capacity (TAWC), used in estimating soil moisture
flow direction, used in the flow accumulation process described below.
Beginning a model iteration¶
The LSM operates on a monthly timestep. At the beginning each step, four state variables must be defined for each pixel:
soil moisture (\(W_s\))
the amount of detained runoff from rainfall (\(D_r\))
the amount of detained runoff from snowmelt (\(D_s\))
the snowpack water equivalent
the number of consecutive months of melting conditions
The three driver variables are used to advance the model for a single onemonth timestep, producing updated versions of the state variables:
These variables are used to perform a water balance on a daily timestep, for each day of the month.
To produce daily water balance inputs from the monthly data, the \(p_{wet}\) parameter is used to divide \(Pr\) among \(n_{wet}\) equallyspaced precipitation days, such that:
Monthly snowmelt is evenly distributed throughout the month.
Daily Water Balance¶
The water balance is an implementation of the Thornthwaite [21][22] equation:
where \(R\) is the rate of surplus water runoff and/or recharge; \(P\) is the effective precipitation rate (equation (2) below); \(E\) is the rate of evapotranspiration (evaporation plus plant transpiration); and \(\frac{dW}{dt}\) is the change in soil moisture.
Effective precipitation (\(P\)) is a function of measured precipitation (\(P_r\)), snow accumulation (\(S_a\)), and snow melt (\(S_m\)):
Note
Snow accumulation and snow melt are represented as snow water equivalent (not snow depth).
Of the quantities in the Thornthwaite’s equation ((1)), only effective precipitation is known directly.
WSIM therefore uses the following steps to arrive at a solution:
Estimate potential evapotranspiration (\(E_0\)), based on day length and temperature.
Use the estimate of \(E_0\), and the known soil moisture, to estimate \(dW/dt\).
Estimate \(E\), based on \(E_0\) and \(dW/dt\).
Solve directly for \(R\).
Potential Evapotranspiration¶
Potential evapotranspiration in the WSIM LSM is calculated using Hamon’s Equation ([8], [9]) to estimate potential evapotranspiration as specified in [24]:
where \(\Lambda\) is the average day length specified as a fraction of the 24hour day between sunrise and sunset, \(T_m\) is the mean temperature in Celsius, and \(e_{T_m}\) is the saturated vapor pressure at \(T_m\) in kPa.
Buck’s Equation ([4][23]) is used to estimate \(e_{T_m}\):
Note
There are numerous formulations for estimating potential evapotranspiration (\(E_0\)) which can, broadly speaking, be divided into two major categories. The first category consists of highly simplified reducedform estimates based on empirical fits for a given reference land cover of short grass. Examples include formulas proposed by Hamon, Thornthwaite, Turc, JensenHaise, Hargreaves, and others ([7][24][16][18][13]). However, it is well known that land cover is a major factor in estimating potential evapotranspiration and it is widely assumed that formulations that take land cover into account are more accurate. All other factors being equal, bare ground will have the lowest potential evapotranspiration and deciduous forest will have the highest.
This gives rise to the second category of formulations that are highly parameterized to include land cover type and many other variables. These processbased or combination methods can explicitly account for different surface characteristics, including vegetation characteristics and the proportion of exposed bare soil. These methods include the PriestlyTaylor, McNaughtonBlack, PenmanMonteith and ShuttleworthWallace methods. The most sophisticated method is the ShuttleworthWallace method, which is a modification of PenmanMonteith (the most commonly used method and the FAO standard). ShuttleworthWallace modifies PenmanMonteith to incorporate a term for bare soil. These methods are all described in [16], [18], and [24]. See also [26] and [27] for a description of ShuttleworthWallace and how it could be parameterized with global data.
Vörösmarty et al. [24] compared 11 different methods of modeling potential evapotranspiration, including methods that either did or did not incorporate differences in land cover. They found that the two best methods for minimizing bias and mean annual error were Hamon’s method and the ShuttleworthWallace method. More recently, Oudin et al. [18] also compared a number of different potential evapotranspiration methods (27 in all). They also found that simple “reference” approaches such as Hamon’s and McGuinness’ performed better than more complex variations. As Oudin et al. [18] wrote: “…if a simple temperaturebased [potential evapotranspiration] estimation works as well as a Penmantype model, why not using [sic] a simpler model with lower data requirements?”
Based on this literature and concurring advice from our science advisors, WSIM chose to implement Hamon’s Equation.
Change in Soil Moisture¶
Change in soil moisture (\(\frac{dW}{dt}\)) is a function of effective precipitation (\(P\) in mm/day), potential evapotranspiration (\(E_0\) in mm/day), soil moisture deficit (\(D_{ws}\) in mm/day), and a unitless soil drying function \(g(W_s, W_c, E_0, P)\).
The soil moisture deficit (\(D_{ws}\)) is the amount of water needed within a time step to fill the remaining soil water holding capacity (\(W_c\) in mm) while satisfying potential evapotranspiration (\(E_0\)). \(W_s\) is the soil moisture in mm.
The unitless drying function, \(g(W_s, W_c, E_0, P)\), is defined as:
The specification follows Vörösmarty et al. [24]. WSIM defines \(g_2(W_s, E_0, P)\) to ensure that when \(P < E_0\), \(g(W_s, W_c, E_0, P) \le W_s\) (i.e., imposing a constraint that \(\frac{dW}{dt} \le W_s\)).
Evapotranspiration¶
Returning to Equation (1), actual evapotranspiration (\(E\)) is calculated as:
Returning to Equation (2), WSIM follows Vörösmarty et al. [24] to estimate snow accumulation (\(S_a\)) and snow melt (\(S_m\)). When monthly average temperature is less than or equal to 1ºC, it assumes all precipitation accumulates as snow pack. This snow pack then melts when monthly average temperature is greater than 1ºC. In elevations less than or equal to 500m, the entire snow pack melts in one month. In elevations above 500m, the snow pack requires two months to melt.
Runoff¶
Finally, WSIM computes two forms of runoff. The runoff as specified above (\(R\)) is always zero during periods when precipitation accumulates as snow pack. This is clearly a falsehood, since most rivers continue to flow in the winter. Therefore, WSIM follows Vörösmarty et al. [24] by including some logic for detention pools (lakes, ponds, shallow groundwater, etc.) that slow down the rate at which runoff as computed above leaves a given grid cell. The revised runoff that accounts for detention pools (\(R'\)) is computed as the sum of detained runoff due to net precipitation (\(R_p'\)) and detained runoff due to snow melt (\(R_s'\)) with a monthly time step as described in Equations (11), (12), and (13) below. (\(D_r\)) and (\(D_s\)) represent the detention pools due to rain and snow, respectively.
where \(z\) is elevation in meters, and \(m\) is the number of consecutive months of melting conditions (\(T > 1 \mathrm{^\circ C}\)).
Flow accumulation¶
After the water balance process has been completed for each day of the month, a flow accumulation algorithm is used to determine the amount of runoff in each grid cell that arrives from upstream locations.
WSIM uses a traditional pixeltopixel based flow accumulation algorithm for this computation. The algorithm uses an eightneighbor flow direction grid that identifies the downstream grid cell for each grid cell. The algorithm has two major benefits: it is well known and easy to implement, and it produces results that highlight the specific paths of the stream network during periods of extreme anomalies.
Flow Direction Specification¶
Flow directions for the pixeltopixel flow accumulator are encoded using the following values:
Direction 
Value 

East 
1 
Southeast 
2 
South 
4 
Southwest 
8 
West 
16 
Northwest 
32 
North 
64 
Northeast 
128 