Abstract

LM3 is a new model of terrestrial water, energy, and carbon, intended for use in global hydrologic analyses and as a component of earth-system and physical-climate models. It is designed to improve upon the performance and to extend the scope of the predecessor Land Dynamics (LaD) and LM3V models by better quantifying the physical controls of climate and biogeochemistry and by relating more directly to components of the global water system that touch human concerns. LM3 includes multilayer representations of temperature, liquid water content, and ice content of both snowpack and macroporous soil–bedrock; topography-based description of saturated area and groundwater discharge; and transport of runoff to the ocean via a global river and lake network. Sensible heat transport by water mass is accounted throughout for a complete energy balance. Carbon and vegetation dynamics and biophysics are represented as in LM3V. In numerical experiments, LM3 avoids some of the limitations of the LaD model and provides qualitatively (though not always quantitatively) reasonable estimates, from a global perspective, of observed spatial and/or temporal variations of vegetation density, albedo, streamflow, water-table depth, permafrost, and lake levels. Amplitude and phase of annual cycle of total water storage are simulated well. Realism of modeled lake levels varies widely. The water table tends to be consistently too shallow in humid regions. Biophysical properties have an artificial stepwise spatial structure, and equilibrium vegetation is sensitive to initial conditions. Explicit resolution of thick (>100 m) unsaturated zones and permafrost is possible, but only at the cost of long (≫300 yr) model spinup times.

1. Introduction

This paper describes a new model, LM3, which is the product of continuing development (Manabe 1969; Milly and Shmakin 2002; Shevliakova et al. 2009) of the land component of the Geophysical Fluid Dynamics Laboratory’s (GFDL) climate and earth-system models. New developments are needed in order to 1) increase the accuracy of means and space–time variabilities of modeled fluxes of water and heat from land to atmosphere and ocean; 2) introduce state variables related more directly to observable quantities, so as to make the model more directly relevant for water-resource impacts and to pave the way for possible future incorporation of water use by humans into the model; and 3) provide spatial detail about subsurface wetness and temperature of relevance to ecological and biogeochemical processes.

Milly and Shmakin (2002) identified flux biases in the Land Dynamics (LaD) model and noted they might be alleviated by vertical resolution of soil-water processes and by explicit inclusion of canopy interception and evaporation therefrom. Subsequent coupling of the model to an atmospheric model suggested that LaD’s complete coupling of ground and canopy temperatures, as well as its neglect of water content in computation of soil thermal properties, led to underestimation of the daily temperature range in the atmosphere near the surface (GFDL Global Atmospheric Model Development Team 2004). Coupling to a physical-climate model (Delworth et al. 2006) revealed errors in seasonal climate of cold regions associated with neglect of phase change of soil water. Swenson and Milly (2006) identified a low bias in the amplitude of annual water storage in the LaD model and tentatively attributed it to the absence of surface water storage in the model. All of these issues are addressed in LM3 presented herein.

The LaD model neglects topographically controlled spatial heterogeneity of soil wetness and hydrologic response, presumably distorting the temporal and spatial variations of runoff, even if parameters are chosen to approximate average fluxes well. Increasingly realistic and physically based approaches have been advanced for addressing such heterogeneity (Wood et al. 1992; Stieglitz et al. 1997; Koster et al. 2000), with use of the TOPMODEL approach (Beven 1986) now being state of the art. TOPMODEL is based on the idea that the water table parallels the land surface and that lateral groundwater flow is concentrated near the surface; these conditions are best met in areas with humid climate and negligible deep permeability. Given this limitation, Koster et al. suggested that TOPMODEL is only a first step in the direction of appropriate treatment of horizontal flows. Herein, we use a generalization of TOPMODEL that extends its applicability to conditions of arid climate and/or high bedrock permeability.

Providing fluxes from land to the atmosphere and oceans is necessary but, increasingly, not sufficient. As the practice of physical-climate modeling evolves, increasing attention is given to analysis of terrestrial impacts of climate change and to interactions between the physical-climate system and the earth’s biogeochemical processes. Model-based investigations in these areas require a land model with expanded treatment of many processes. Representation of biophysical and biogeochemical processes requires detailed physical information, for example, on the vertical distribution of water, ice, and temperature within the soil profile. Additionally, the water-resource relevance of a model is enhanced by explicit treatment of groundwater and surface water stores, including lakes and rivers.

Given these considerations, we developed a new land model, LM3, which includes enhanced conceptualizations of snowpack, subsurface water, canopy biophysics, rivers, and lakes. To support more realistic treatment of canopy biophysics, we incorporated the dynamic vegetation model LM3V of Shevliakova et al. (2009).

Objectives of this paper are to outline the structure of LM3 and to describe its behavior by comparison with observed data. We describe and evaluate two model configurations—LM3.0, used in phase 5 of the Coupled Model Intercomparison Experiment (CMIP5), and a refined version, LM3.1.

2. Model description

a. Overview

In LM3, global land area is divided into a grid of cells; cell boundaries may line up with latitude and longitude lines or be based on any other local row–column coordinate system, as in a cubed-sphere grid. Each cell area is decomposed into one or more tiles, within which resolved processes are treated as horizontally homogeneous, that is, varying only in the vertical direction. Each tile has a substrate of soil, lake, or glacier. Total area of any substrate within a cell does not change. Any lake or glacier is permanently represented by a single tile, but the number of soil tiles is not necessarily constant. Soil tiles (only) can contain vegetation. Vegetation disturbance (e.g., harvest and fire) can be applied to an area smaller than the total vegetated area in a cell, and this is accommodated by dynamic tiling of soil, with a new tile introduced for each new vegetation cohort. Soil tiles are merged when growth causes properties of previously disturbed vegetation to converge. Tiles can also be used to represent other types of subgrid variability, for example, with respect to soil properties or topographic relief, but this feature is not exercised in this paper; rather, it is assumed that hillslope-scale topography is the major control of spatial variation of soil hydrologic response.

Each tile interacts with the atmosphere, independently of other tiles, through the atmospheric surface boundary layer. Atmospheric input to LM3 consists of precipitation; downward short- [visible (VIS) and near infrared (NIR), direct and diffuse] and longwave radiation; and air temperature, humidity, pressure, and wind speed. When the land model is coupled to an atmospheric model, the resultant surface fluxes are fed back to the atmospheric model. A transient snowpack may be present in any tile, regardless of substrate. Where present, the snowpack (and not the substrate) is the site of all ground-level energy exchange between the tile and the overlying air. Deep snowpack can mask a fraction of the plant canopy; fraction visible to sky decays with increase of snow depth.

Soil and glacier tiles can produce runoff. Runoff is deposited into a river network whose reaches are resolved on the land grid. Rivers are conceptualized as having zero surface area and, hence, do not interact directly with the atmosphere. However, rivers are coupled horizontally with lake tiles, which have their own water and energy balances, so indirect atmospheric interaction with rivers does occur.

The sensible heat capacity of water and ice cannot be ignored in lakes, snow, or glaciers. Thus, in order to conserve energy, the sensible heat capacity of water is included in all stores and fluxes throughout the model.

For all three types of substrate, vertical mass flux across the bottom boundary is specified as zero, and heat flux is set to the geothermal heat flux. (Numerical values of this flux and other model inputs are configuration dependent and are given in section 3.) The top surface of all three substrates is opaque, and reflectance is characterized with a bidirectional reflectance distribution function (BRDF; Schaaf et al. 2002) for VIS and NIR wave bands.

b. Vegetation

Carbon and vegetation dynamics are modeled as described by Shevliakova et al. (2009). For each tile, the vegetation model dynamically computes 1) a single vegetation type (C3 grass, C4 grass, temperate deciduous tree, tropical tree, or cold evergreen tree) as a function of climate and vegetation biomass (where the latter generally varies across tiles in the same cell, as a result of disturbance) and 2) physical structure of above- and belowground biomass, including leaf and root phenology, on the basis of allometric equations. Carbon is tracked in five stores: leaves, fine roots, heartwood, sapwood, and labile. Carbon dynamics are driven by weather and climate, atmospheric CO2 concentration, soil state, and natural and anthropogenic disturbances.

Net photosynthesis is determined mechanistically (Farquhar et al. 1980; Collatz et al. 1991, 1992). Stomatal conductance is set to the smaller of a non-water-stressed value and a water-limited value. The non-water-stressed value is parameterized following Leuning (1995), but with cuticular conductance neglected. By making conductance directly proportional to leaf area, this approximation allows an analytical solution of the coupled equations with consideration of change in photosynthetic limitation (i.e., by light, RuBisCo, carboxylation, or CO2) vertically through the canopy. The water-limited conductance is that which balances transpiration with the hypothetical rate of water uptake from soil–root interface to leaves whose water potential is at the wilting point. This potential-driven flow is determined by plant structural parameters (root density, plant height, and sapwood mass), root permeability , and xylem resistance . Negative uptake (and associated “hydraulic lift”) can optionally be permitted. Leaf drop can be triggered by drought or cold stress.

For radiative transfer, the canopy is treated as a dispersed medium, and the two-stream approximation is used; leaf reflectance is adjusted for snow interception. Canopy interception and heat storage are treated with a single lumped store (Deardorff 1978) whose water capacity is proportional to leaf area index (LAI; Dickinson et al. 1993). The canopy air space, which is characterized by lumped stores of temperature, humidity, and CO2 concentration, interacts through aerodynamic resistances with the canopy, the atmosphere above, and the soil or snow below.

c. Snowpack

Bulk density and thermal conductivity of snow are constant; metamorphosis is ignored. Pack grows by accumulation of snowfall and shrinks by surface sublimation and melting/drainage. Snow storage in excess of 1000 kg m−2 is converted to frozen runoff to the river system. The pack is divided into a constant number of layers that all expand and contract with the changing snow mass. Each layer has a capacity to hold liquid water; excess liquid flows out the bottom of the layer. Ventilation, internal vapor transport, and metamorphosis are ignored.

The shortwave reflectivity of the snowpack varies continuously from that of snow substance (assumed constant; aerosol deposition is ignored) when the snow is much deeper than some masking depth (a model parameter) to that of bare soil as the pack depth approaches zero.

d. Soil and bedrock

Soil and bedrock are conceptualized together as a porous medium having depth-dependent physical properties. The vertical domain is discretized into an arbitrary number of layers whose thickness generally increases with depth. Water movement is treated with a dual-domain model (Beven and Germann 1982); in parallel with (but not in equilibrium with) a ubiquitous microporous (Darcian) domain, a macropore network is present wherever heartwood biomass (a stable measure of biological activity) exceeds some threshold . The macropores instantaneously transport would-be infiltration-excess runoff (see below) into the soil, filling any available pore space from a specified depth upward; storage in macropores is ignored.

In the Darcian domain, liquid water conservation is expressed in one-dimensional (z positive downward) form, with lateral flow divergence G, uptake by roots r, and absorption of water from macropores M parameterized as vertically distributed sinks/source:

 
formula

with

 
formula

where t is time, is density of liquid water, is volumetric liquid water content, is mass flux of water, is rate of freezing, K is vertical hydraulic conductivity, and is matric head of water (Hillel 1980). Energy is stored as sensible and latent fusion heat and is transported by conduction and flowing liquid:

 
formula

with

 
formula

where is density of ice, is volumetric ice content, T is soil temperature, C is bulk soil heat capacity, is latent heat of fusion of water, is heat flux, is soil thermal conductivity, is specific heat of water, and TM is temperature of water absorbed by macropores.

At the bottom of the soil–bedrock column, a zero water flux and a constant geothermal heat flux are prescribed. At the surface, water flux is tentatively set to the difference between any liquid input reaching the ground from rainfall, snowmelt, or the interception store, minus evaporation. However, if the resultant solution would cause matric head at the surface to become positive (infiltration capacity of soil exceeded by water supply rate), the matric head is instead set to zero (thinly ponded surface condition), the continuity equation is again solved, and water that does not infiltrate is tentatively defined as infiltration-excess runoff, which will, however, enter the soil if macropores are present and pore space is available. The heat flux at the soil surface is determined by balance of the radiative and turbulent fluxes, with consideration of latent heat of evaporation/sublimation and sensible heat advection by water fluxes. Turbulent fluxes between the surface and the overlying air are determined with consideration of effects of vegetation on the aerodynamic resistances. Water vapor flux from the soil surface to the overlying air is impeded additionally by a surface resistance associated conceptually with the presence of a litter layer, but litter is not explicitly modeled. Rather, the litter resistance to vapor transport is made directly proportional to the heartwood biomass, which is used as a surrogate for litter mass. Relative humidity h of air at the soil surface (below any litter) is determined by the thermodynamic relation (Edlefsen and Anderson 1943, p. 145):

 
formula

in which g is the acceleration of gravity and is the gas constant for water vapor.

At the soil surface, under unsaturated, unfrozen conditions, nonhysteretic power laws relate , , and K (Cosby et al. 1984; Campbell 1974):

 
formula

and

 
formula

where θ0, ψ0, B, and K0 are soil constants. The soil-freezing curve is approximated as a step function: all phase change occurs at a soil-water-freezing point, which is a soil constant. Where liquid and solid water coexist, matric potential of liquid determined for unfrozen conditions is scaled by a surface-tension ratio (liquid against ice versus liquid against air) of 2.2 (R. D. Miller 1980). The difference in density of water between liquid and frozen states is ignored, and the soil–bedrock is considered rigid and unaffected by freeze–thaw processes.

The vertical variation of soil physical properties is treated by assuming microscopic geometric similarity, viscous flow, and dominance of matric head by surface tension (E. E. Miller 1980). It follows that the relations among , , and K at any depth can be determined from the relations at a single reference depth and a depth-dependent microscopic length factor , and that the dependences of and C on soil-water content are independent of depth. Specifically,

 
formula

and

 
formula

The microscopic length scale varies with depth according to

 
formula

which means that the value of K at saturation decays exponentially with depth. The e-folding depth conceptually represents the soil thickness, and is the ratio of saturated hydraulic conductivity at large depth (i.e., in bedrock) to that at the surface .

To determine r, roots are represented as cylinders of small radius, and difference between bulk water potential and water potential at the soil–root interface is determined by the near-field steady-state solution of the flow equation (Gardner 1960).

To solve (1), the lateral flow divergence must also be specified. Note that the vertical integral of G, which we denote by , is the rate of subsurface water discharge to surface water. To set the stage for evaluation first of , we conceptualize a tile as having a characteristic hillslope structure. We let x denote the distance from stream to any location on the hill, with l being the distance to the top of the slope, and we describe the shape of the hill with an elevation (above stream level) function Z(x) and a width function w(x), the latter being the width of the hill per unit length of stream. Let D(x) be the depth of the water table below the surface. We consider a steady-state groundwater (i.e., below the water table) flow associated with a uniform groundwater recharge (equals discharge at steady state) rate . The horizontal groundwater flow at any x, equal to the accumulation of recharge uphill from x, can be expressed, according to Darcy’s law, as the product of the transmissivity R (vertical integral of horizontal hydraulic conductivity below the water table, a function of water-table depth), the slope width, and the horizontal hydraulic gradient, which is assumed (under the “Dupuit approximation”) to be given by the slope of the water table. The resultant mass-balance equation is

 
formula

However, where and when (11) would imply a water table above the surface, D is set to zero, implying presence of a seepage face. From (10), and ignoring anisotropy of hydraulic conductivity, we would have

 
formula

in which is an effective bedrock thickness potentially contributing to horizontal hillslope-scale flow. For the groundwater solution, we approximate and augment this relation to obtain

 
formula

in which b is effective bedrock thickness and is enhancement of lateral saturated hydraulic conductivity associated with macroporosity and anisotropy. (However, whereas macropores permit infiltration into frozen soil, it is assumed they do not permit relatively long-distance lateral transport of that infiltrated water through frozen ground.) The combination of (11) and (13) can be viewed as a generalization of TOPMODEL (Beven 1986), to which it reduces, in essence, when and horizontal gradients of D are neglected (e.g., for relatively impermeable bedrock, steep slopes, and a humid climate).

The continuity Eq. (11) is integrated numerically in dimensionless form for the possible ranges of three relevant parameters: , , and , in which is , as well as the given shape functions and . The solution gives dimensionless water-table depth as a function of , which can be averaged over the slope to determine average water-table depth. The solution also identifies what areal fraction of the slope has zero water-table depth. Thus, we can construct lookup tables that define recharge and saturated-area fraction as functions of average water-table depth. When the vertical flow in Eq. (1) is to be solved, the water-table depth is first determined from the profile of , and then groundwater discharge and saturated area are found from the lookup tables; if the water table is absent from the modeled domain, its depth is estimated by assuming hydrostatic conditions between the domain bottom and the water table. Effective rainfall onto the saturated-area fraction is rejected from infiltration into the soil column and becomes surface runoff. The groundwater discharge is distributed vertically into G in proportion to the contributions of layers to transmissivity. Where permafrost is present, however, is assumed to be zero; at such locations, lateral flow from the active layer is evaluated as the product of average elevation gradient and unsaturated hydraulic conductivity.

e. Lakes

Lakes are characterized by changing depth d and vertical profiles of temperature and ice content. Lakes (including lake ice) receive water from melting of snowpack above lake ice, from precipitation, and from river inflows. For large lakes that span more than one grid cell, lake levels are kept equal across cells. Lakes lose water by evaporation and by outflow to rivers Q given by the stage–discharge relation for an ideal broad-crested weir (Henderson 1966):

 
formula

in which W is an effective width of lake at outflow and d0 is the value of d at which outflow ceases. Lake mass changes are applied at the surface. Lake bottoms are impermeable but receive the geothermal heat flux. Lake temperatures change as a result of stability-dependent vertical diffusion, following Henderson-Sellers (1985) and Hostetler and Bartlein (1990), but with addition of a constant to the vertical eddy diffusivity to crudely parameterize the unresolved effects of three-dimensional mixing in large, deep lakes. Ice and liquid are moved vertically as needed to keep ice above liquid. When lake level falls below the level of outflow, evaporation is artificially reduced from its wet-surface value by a factor :

 
formula

in which p is the value of d/d0 at which evaporation ceases. This scaling serves as a surrogate for reduction of lake surface area. Lake storage is prevented from growing unbounded by requiring that all lakes have an outlet; thus, all land areas potentially drain to the oceans. Over time, lake layers are split or merged as necessary to maintain a constant number of layers whose thicknesses do not differ by more than a factor of 2.

f. Rivers

Every land cell has a single river reach, which transports water (liquid or solid) into a downstream cell or into the ocean. Any runoff produced by soil and/or glacier tiles in the cell is combined with water flowing into the cell from any upstream river reaches and is first fed into a lake tile, if one is present in the cell. Outflow from the lake enters the cell’s river reach. This approach is even used in multicell lakes, in which case the river network prescribes a pathway that snakes through all relevant cells and the river reaches are artificial. Rivers do not lose water to soil tiles; neglected are such processes as streamflow losses to infiltration in arid lands, transient bank storage, and infiltration of overbank floodwater into soil.

River water changes phase as appropriate at the freezing point, but, for simplicity, ice flows as freely as liquid. Except where backwater (downstream control of flow) is present, the “at-a-site” and “downstream” hydraulic geometric relations of Leopold and Maddock (1953) together yield the reach discharge as a function of reach storage S and long-term-mean discharge :

 
formula

in which and MR are the downstream coefficient and exponent for velocity, m is the at-a-site exponent for velocity, and is the length of the river reach. The long-term mean, which represents channel-forming processes, is evaluated by a 10-yr exponential smoother. Presence of backwater is prescribed in any reach that is separated from the ocean by one or more other reaches, all of which have a topographic slope lower than some global parameter. In backwater reaches, river and lake levels rise and fall with the water level of the next reach downstream.

g. Glacier

Glacier refers to the Greenland and Antarctic ice sheets and to mountain glaciers that may be resolved by the model. Glacier presence is specified as input to the model. Treatment of glaciers follows the very simple approach used in the LaD model. Glaciers store and conduct sensible heat but are not allowed to change mass; neither runoff nor sublimation is treated. Melt is allowed, but meltwater is permanently sequestered below the surface of the glacier, never to refreeze. However, snowpack that lies on top of glaciers interacts with the atmosphere, as already described. Runoff from the snowpack is sent through the river network in the same manner as liquid runoff.

h. Overview of solution procedure

For each time step, the vertical radiative, diffusive, and conductive fluxes and the substrate–surface melt are computed simultaneously and implicitly, with coupling through any snowpack and through the whole substrate column. Next, movement of water in the condensed phases is computed in downstream fashion from precipitation to canopy to snowpack to soil (or lake or glacier) to rivers to ocean discharge, with adjustment of temperatures for associated advection of heat. When this computational chain reaches the soil column, the distributed sinks for groundwater divergence and uptake by plants are evaluated explicitly and then the solution for flow in the full column is determined implicitly.

3. Experiments

Numerical experiments to evaluate LM3 performance were run with seven combinations of model configurations (LM3.0 and LM3.1), initial conditions (base, DRY, and GREEN), and atmospheric boundary conditions (base and snow S), all described below. Combinations used were LM3.0, LM3.0-DRY, LM3.0-GREEN, LM3.1, LM3.1-DRY, LM3.1-GREEN, and LM3.1-S.

a. Boundary and initial conditions

The LM3.0 and LM3.1 experiments consisted of a spinup phase followed by a historical phase. In the 270-yr spinup phase, atmospheric forcing was constructed by cycling 9 times through the historical data (Sheffield et al. 2006) for 1948–77, atmospheric CO2 concentration was set to the preindustrial level (286 ppm), and vegetation disturbance by land use was neglected. In the 331-yr historical (1678–2008) phase, which was initialized from the end of the spinup phase, atmospheric forcing cycled 10 times though the 1948–77 data, and then the 1978–2008 data were used; historical atmospheric CO2 concentrations (Meinshausen et al. 2011) and vegetation disturbance by land use (Hurtt et al. 2011) were prescribed. Geothermal heat flux was zero in LM3.0 and 0.065 W m−2 in LM3.1 (Pollack et al. 1993). Subsurface water was initially in hydrostatic equilibrium, with a water-table depth of 10 m and a temperature of 288 K. Rivers were dry, snow was absent, and lakes had d at d0 and temperature at 288 K. Glacier temperature was 260 K. Vegetation biomass was zero.

Sensitivities of the LM3.0 and LM3.1 experiments to initial conditions were explored to assess system time scales and any dependence of spun-up state and fluxes on initial conditions. In the LM3.0-DRY and LM3.1-DRY experiments, the initial water-table depth was changed from 10 to 100 m. In the LM3.0-GREEN and LM3.1-GREEN experiments, the initial biomass stores were changed from zero to globally constant values more than double the largest equilibrium values generated in LM3.0 and LM3.1 experiments. The DRY and GREEN experiments were run only for the 270-yr spinup phase.

The 61-yr (1948–2008) LM3.1-S experiment, intended to assess sensitivity to uncertainty in frozen precipitation, was branched from the historical phase of the LM3.1 experiment at the start of year 1948; it differed from LM3.1 only in that all snowfall was increased by 50%. An adjustment of this magnitude was motivated by the spatially biased (low elevation) sampling of snowfall in remote areas of high relief and was found to be supported, for example, by the precipitation dataset used by Milly and Shmakin (2002). As we shall see, some adjustment of precipitation amounts appears necessary in some northern basins to reconcile precipitation observations with measurements of streamflow.

b. Discretization

All experiments used the 2° latitude × 2.5° longitude grid that was used for the Earth System Model with the Modular Ocean Model (ESM2M; Dunne et al. 2012) in its application for CMIP5. Most processes were performed on a 30-min time step; river and horizontal lake calculations were performed on a 1-day time step. LM3.0 soil–bedrock domain depth was 10 m, with 20 layers ranging from 0.02 to 2.5 m in thickness. LM3.1 depth was 200 m, with 35 layers ranging from 0.02 to 20 m. Glacier domain was 6 m, in 18 layers ranging from 0.02 m at the surface to 1 m at the bottom. Lakes were divided into 20 dynamic layers. Snow was divided into five dynamic layers accounting sequentially for 5%, 20%, 50%, 20%, and 5% of the pack.

c. Model configurations

The two model configurations (LM3.0 and LM3.1) differ in many respects (Table 1). The LM3.0 configuration is the same as that used for integrations with the GFDL ESM2M (Dunne et al. 2012); simple variations on LM3.0 were also used in other GFDL models contributing to CMIP5. The parameter settings for LM3.0, which are not exhaustively documented here, were based mainly on a combination of values from LM3V (Shevliakova et al. 2009) and the LaD model (Milly and Shmakin 2002) and, for features new to LM3, on literature estimates and subjective judgment.

Table 1.

Parameter settings for the two model configurations.

Parameter settings for the two model configurations.
Parameter settings for the two model configurations.

The LM3.1 configuration is the result of model evolution subsequent to the freezing of LM3.0 for CMIP5. Some changes were motivated by deficiencies in LM3.0 performance; in such cases, the motivation is noted herein. Others were made in order to refine relatively ad hoc parameter values on the basis of literature search and analysis subsequent to freezing of LM3.0, whether or not they had significant impact on model performance.

For both model configurations, tiling of glacier, lake, and soil was determined by a merger of data from Zobler (1986) with lake areas defined as “water body” in the U.S. Geological Survey (USGS) Global Land Cover Characteristics database (http://edc2.usgs.gov/glcc/glcc.php). Grid cells were not tiled with respect to differing soil types or slopes. A drainage route to the ocean was defined for all land cells. The river network was regridded from the 30-min University of New Hampshire Simulated Topological Network (STN-30) to the model grid by means of automated procedures (see, e.g., Fekete et al. 2001) and further modified manually with reference to Barraclough (1988).

For LM3.1, the BRDF parameters were estimated a priori (i.e., through a process that did not use the land model) from the Moderate Resolution Imaging Spectroradiometer (MODIS) global BRDF parameter data product (MOD43C2, version 4). For soil, the parameters were assigned a minimum-variance weighted average of two values, with weights depending on annual-mean precipitation (a surrogate for vegetation density); the two values, both derived under snow-free conditions, are a global constant value (used where vegetation is dense and soil is not visible to MODIS) and the local observed MODIS value, the latter adjusted approximately for partial cover by vegetation. For other substrates and snowpack, the parameters were assigned globally constant values estimated subjectively from visual displays of MODIS values over representative target areas.

For LM3.1, leaf reflectances were set for the tropical and cold evergreen trees so that areal median differences between model and MODIS (global albedo data product MOD43C1, version 4) estimates of time-average surface reflectance for diffuse VIS and NIR radiation were zero. Data were masked to include only areas and months with air temperature above 285 K (to exclude snow effects) and LAI greater than 2 (to minimize soil influence). Leaf-reflectance parameters for grasses and for deciduous trees were then set equal to those for the tropical trees. Subsequently, snow-masking depths for cold evergreen forest were estimated in a similar manner, but using cold-season data. For other vegetation types, the masking depth was made very large; its value is irrelevant both in the leaf-off season and in tropical forests, where snow does not occur.

In LM3.1, the liquid interception capacity per unit leaf area was set on the basis of field observations in the range 0.1–0.2 mm (Brutsaert 2005, p. 106), with the higher value pertaining to (mostly needleleaf) cold evergreen forests (Bonan 2002, p. 37) and the lower value to all other vegetation types. These values were tripled for interception of snow (Bonan 2002, p. 136).

When LM3.0 was formulated, xylem resistance was neglected. Under an ecological optimality hypothesis (Milly and Dunne 1994), root permeability was subjectively set just high enough that modeled global evapotranspiration was minimally sensitive to the permeability. In LM3.1, root resistance was assumed negligible (large permeability), and its effect was instead lumped into the xylem resistance. The globally constant xylem resistance was adjusted approximately to minimize errors in modeled annual runoff from a set of 53 river basins. The value so determined was within the central range of measured values of xylem resistance (Maherali et al. 2004).

Drought-induced leaf drop in LM3.0 follows the rule used in LM3V: leaves drop when monthly average , exponentially weighted in depth (with an e-folding depth on the order of 1 m), falls below a value prescribed as a function of vegetation type (0.15 for C4 grass and deciduous forest and 0.3 for other types) and independent of soil type. In LM3.1, according to a more physically based approach, drop is triggered when monthly average , vertically weighted by root biomass, exceeds a permanent wilting point of 150 m. Cold-triggering of leaf drop (except in cold evergreen trees) in LM3.0 follows LM3V: drop occurs at 5°C for C3 grass and 10°C for other (except cold evergreen) vegetation types. In LM3.1, the critical temperature is 0°C. This change was made because it was found to improve the modeled seasonal march of green-up, as compared to MODIS reflectances, in some midlatitude regions.

For all but the cold evergreen–temperate deciduous forest boundary, vegetation boundaries in both configurations are determined dynamically according to biogeographical rules defined by Shevliakova et al. (2009) for LM3V. The cold evergreen tree type is specified when the number of cold months (mean air temperature below 10°C) is in the range 9–12 for LM3.0 or 7–9 for LM3.1. The change from LM3.0 (LM3V) to LM3.1 was made in order to better locate the region of cold evergreen forest, as reflected in MODIS reflectance observations.

In LM3.0, snow density and thermal conductivity were assigned geometric means of values for “old” and “new” snow given by Garratt (1994, p. 291). This combination was later found to be physically inconsistent, with excessive thermal conductivity (Sturm et al. 1997), and was replaced in LM3.1 by a consistent combination based on a density of 250 kg m−3, which is the central value for “settled snow” given by Cuffey and Paterson (2010).

For both model configurations, , , , and B were assigned texture-specific values based on Zobler (1986) and Cosby et al. (1984). Thermal properties were assigned on the basis of measurements given by Brutsaert (1982) and Garratt (1994). The soil-water freezing temperature was given a globally constant value of 2 K below the freezing point of pure water; this is near the middle of the range of values typical for freezing of loosely bound water in soils (Lunardini 1981, p. 103).

For simplicity, hill slopes were treated as rectangular ramps in both models. Thus, w(x) = 1 and Z(x) = . Hillslope length l was set to 1000 m everywhere. For LM3.0, was assigned values of 40, 160, or 320 m, according to slope classes of flat, medium, or steep in the Digital Soil Map of the World and Derived Soil Properties (FAO 2003) and then spatially averaged to the model grid. These values were halved for LM3.1, on the basis of typical relations between slope gradient and relief intensity (FAO 2006).

Runoff dynamics are controlled jointly by a poorly constrained set of parameters, including soil thickness scale, effective bedrock thickness, macropore conductivity, and bedrock hydraulic conductivity. These were set subjectively by a combination of 1) a priori assignments and 2) sensitivity analysis of temporal variability of daily river discharges and spatial extent of saturated areas. An upper limit on b for deep, isotropic bedrock can be deduced from the solution of Tóth (1962) to be l/π, or about 300 m; values for finite bedrock and for media with relatively small vertical permeability could be much smaller. In LM3.0, where b, , and were zero by definition, it was necessary to give a large value (40 m) in order to have significant lateral subsurface flow; otherwise, soils became waterlogged, and saturated areas grew very large to produce all runoff as saturation-excess runoff. For LM3.1, effective bedrock thickness and soil thickness scale were assigned globally constant values of 10 and 0.5 m, respectively; Kb was computed from the permeability map of Gleeson et al. (2011); and macropore hydraulic conductivity was set to 10−3 m s−1. With larger values of zs and/or b, the modeled seasonal cycle of discharge was too weak. With much smaller macropore conductivity, soils waterlogged and saturated areas expanded to unrealistic proportions. Threshold wood biomass for presence of macropores in LM3.1 was set so that macropores generally would be present in humid regions and absent in arid regions.

The litter resistance coefficient was zero in LM3.0. With no resistance, evaporation from soil was an unrealistically large fraction (15%–30%) of total evapotranspiration from forests (see, e.g., Stoy et al. 2006). In LM3.1, the coefficient was given a globally constant value that resulted in evaporation from soil being a small fraction (5%–10%) of total evapotranspiration from forests.

The critical river slope for backwater was zero in LM3.0 and was chosen to produce backwater on most of the Amazon River in LM3.1; this choice resulted in extensive backwater also in the Parana and Ob basins. For the 16 large lakes, outflow widths were arbitrarily set to 200 m in LM3.0 and 800 m in LM3.1, and lake depths were assigned as lake volume divided by surface area, according to data from van der Leeden et al. (1990), but no greater than 50 m in LM3.0. In all other lakes, transient storage above d = d0 was prohibited, lake area fraction was defined on the basis of the USGS Global Land Cover Characteristics database “IGBP Water Bodies” field, and d0 was set to 2 m. A background lake eddy diffusivity was introduced for LM3.1 and set by optimization with observations. Lake surface roughness length changed from LM3.0 to LM3.1 after an unintended value of 0.01 m was found to have been assigned in LM3.0. For LM3.0, parameters describing hydraulic geometry of rivers were assigned central values from literature covering more than 20 multisite field studies (e.g., Leopold et al. 1964; Emmett 1972; Ibbitt 1997; Wohl 2004) and a survey of data from more than 200 sites (Xu 2004). For LM3.1, parameters were set close to those from the original work of Leopold and Maddock (1953), because this resulted in more realistic simulation of the annual range of total water storage.

4. Results and analysis

a. Spinup and mass and energy balances

In LM3.1, which has a 200-m vertical domain, global averages of total water storage, frozen water storage, and temperature are still changing after 270 years of spinup (Fig. 1). The difference in water storage between LM3.1 and LM3.1-DRY, associated mainly with water-table depth, has disappeared in humid regions (Fig. 2), where initially low LM3.1-DRY water tables rise rapidly by groundwater recharge from above. In semiarid and arid regions, however, the fall of the initially high water table in LM3.1 is limited by slow rates of groundwater discharge to streams, and the rise of initially low water table in LM3.1-DRY is limited by small recharge from the surface. Linear extrapolation of the difference in global-mean rates of storage at 270 years implies time to water-storage convergence, hence equilibrium, on the order of 2000 years. However, initial differences in near-surface water content will literally be frozen permanently in cold regions.

Fig. 1.

Annual spinup time series of averages over all land area of water storage, frozen water storage, depth-average soil–bedrock temperature, and biomass for differing initial conditions of LM3.1.

Fig. 1.

Annual spinup time series of averages over all land area of water storage, frozen water storage, depth-average soil–bedrock temperature, and biomass for differing initial conditions of LM3.1.

Fig. 2.

Maps of water-table depth (m) after 270-yr spinup of (top) LM3.1 and (bottom) LM3.1-DRY. Contour lines are at 3, 10, and 30 m.

Fig. 2.

Maps of water-table depth (m) after 270-yr spinup of (top) LM3.1 and (bottom) LM3.1-DRY. Contour lines are at 3, 10, and 30 m.

Changes in temperature and frozen water storage at depth are limited by heat conduction, and the character of conduction is apparent in the LM3.1 ice time series (Fig. 1) because soils were unfrozen everywhere initially. The global-mean temperature curves do not go as , reflecting systematic differences in thermal diffusivity between warming and cooling latitudes, related to differences in water and ice content of the conducting medium. Simple theory suggests time to thermal equilibrium for a 200-m column is thousands of years. After the 270-yr spinup, however, a near–steady state, with a realistic geothermal temperature gradient of ~0.02 K m−1, has been established at 40°N, where the initial temperature was not far from the equilibrium value (Fig. 3).

Fig. 3.

Depth–latitude dependence of temperature during LM3.1 spinup. Blue line indicates soil-freezing temperature.

Fig. 3.

Depth–latitude dependence of temperature during LM3.1 spinup. Blue line indicates soil-freezing temperature.

Global biomass in LM3.1 stabilizes relatively quickly (Fig. 1). After about 200 years, however, the global difference between LM3.1 and LM3.1-GREEN ceases its decrease. The residual difference is associated with persistent differences in locations of boundaries between vegetation types (Fig. 4). In LM3.1-GREEN, the deciduous and tropical forests extend into areas (total area 1.4 × 107 km2) that are C3 and C4 grasslands, respectively, in LM3.1. The impacts on biophysical properties and water and energy balances are substantial.

Fig. 4.

Maps of vegetation type after 270-yr spinups of (top) LM3.1 and (bottom) LM3.1-GREEN. Hatching in LM3.1-GREEN map indicates areas where vegetation type is the same as in LM3.1; solid coloring indicates where vegetation type differs.

Fig. 4.

Maps of vegetation type after 270-yr spinups of (top) LM3.1 and (bottom) LM3.1-GREEN. Hatching in LM3.1-GREEN map indicates areas where vegetation type is the same as in LM3.1; solid coloring indicates where vegetation type differs.

In spite of the sensitivity of model states to initial conditions on water, the transient rates of water and energy storage in LM3.1 after 270 years are small in comparison with the surface fluxes (Table 2), and the global fluxes are relatively insensitive to the differences in initial conditions. At 270 years, global-mean evapotranspiration (latent heat flux) and runoff both differ across initial conditions by only about 2 mm yr−1 (0.15 W m−2). Thus, the land can be considered spun-up from the perspective of simulation of physical climate. Note, however, that the LM3.1 water-storage rate at 270 yr (−1 mm yr−1) is still sufficient to drive a spinup-related sea level rise on the order of 0.5 mm yr−1. The negligible mass- and energy-balance residuals confirm the conservation of mass and energy by the model.

Table 2.

Water and energy balances of global landmass, expressed per unit area of land, for years 241–270 of LM3.1 spinup. (Sums and residuals were computed before rounding addends.)

Water and energy balances of global landmass, expressed per unit area of land, for years 241–270 of LM3.1 spinup. (Sums and residuals were computed before rounding addends.)
Water and energy balances of global landmass, expressed per unit area of land, for years 241–270 of LM3.1 spinup. (Sums and residuals were computed before rounding addends.)

The relatively shallow (10 m) soil–bedrock domain in LM3.0 allows more rapid spinup (not shown) than in LM3.1. After 270 years of spinup, the water-storage rate in LM3.0 is approximately −0.1 mm yr−1, and temperatures have reached equilibrium throughout the column. Biomass in LM3.0 spins up at the same rate as in LM3.1, and the sensitivity of equilibrium biomass to initial conditions seen in LM3.1 is also seen in LM3.0

b. Vegetation

Modeled vegetation type is sensitive to model configuration. In LM3.0 (Fig. 5), relative to LM3.1 (Fig. 4), the cold evergreen forests are displaced northward to the Arctic Ocean coast and replaced in the south by temperate deciduous forests and C3 grasslands; these shifts are a consequence of the different values of biogeographical parameters. Because LM3 biomes correspond to only five single-plant functional types, we do not attempt a comparison with observed biomes, which have much greater variety. Instead, modeled and observed vegetation are compared by use of a normalized difference vegetation index (NDVI). As an overall measure of vegetation density, NDVI is dependent upon climate, soil, natural vegetation dynamics, and land use. Additionally, NDVI is strongly affected by snow cover; we focus here on snow-free conditions.

Fig. 5.

Map of vegetation type after 270-yr LM3.0 spinup. Hatching indicates areas where vegetation type is the same as in LM3.1; solid coloring indicates where vegetation type differs.

Fig. 5.

Map of vegetation type after 270-yr LM3.0 spinup. Hatching indicates areas where vegetation type is the same as in LM3.1; solid coloring indicates where vegetation type differs.

The broadscale spatial structure of LM3 vegetation density during northern summer in both model configurations is similar to the observed, but systematic discrepancies are notable (Fig. 6). In contrast to MODIS, both configurations of the model produce NDVI that is relatively constant within biomes, with relatively narrow transition regions between. LM3 does not reproduce the observed peaking of NDVI values above 0.8 in the centers of the cold evergreen, temperate deciduous, or tropical forests, but rather produces broad “plateaus” of vegetation density. Conversely, NDVI tends to be too high outside the forests, specifically in the subtropics and the Tibetan Plateau. The greater positive NDVI bias of LM3.0 in far northern Eurasia and North America is associated with the presence there of cold evergreen forest in LM3.0. Differences between LM3.0 and LM3.1 in Australia and southern Africa are associated with differences in parameterization of leaf drop.

Fig. 6.

Maps of June–August average of NDVI, for MODIS, LM3.0, and LM3.1, along with model–MODIS difference maps. NDVI is computed for each month as the difference between surface NIR and VIS reflectances, divided by their sum. NIR and VIS reflectances are averaged over 2001–06 for MODIS and 1999–2008 for the model. Here and in other figures, mn and sd are area-weighted mean and std dev of observations; mnΔ and rmsΔ are area-weighted mean and RMS values of difference between model and observations; and r is correlation between model and observations.

Fig. 6.

Maps of June–August average of NDVI, for MODIS, LM3.0, and LM3.1, along with model–MODIS difference maps. NDVI is computed for each month as the difference between surface NIR and VIS reflectances, divided by their sum. NIR and VIS reflectances are averaged over 2001–06 for MODIS and 1999–2008 for the model. Here and in other figures, mn and sd are area-weighted mean and std dev of observations; mnΔ and rmsΔ are area-weighted mean and RMS values of difference between model and observations; and r is correlation between model and observations.

c. Surface albedo

Because snow is highly reflective, biases in albedo are most substantial in cold regions (Fig. 7). Effective annual albedo in cold regions is dominated by that in springtime, when snow is present and solar radiation is substantial, so we focus on northern spring. LM3.0 has a prominent zonal pattern of bias in cold regions—a broad band of high bias across the high middle latitudes, with a band of low bias at even higher latitudes. This error pattern is absent in LM3.1, which nevertheless displays substantial regional biases in cold regions. The improvement in albedo from LM3.0 to LM3.1 is a direct consequence of the differing locations of cold evergreen forests (Fig. 5).

Fig. 7.

Maps of March–May average of effective white-sky (diffuse) land albedo, model (1999–2008) minus MODIS (2001–06), for (top) LM3.0 and (bottom) LM3.1. Effective albedo is ratio of time averages of up- and downward shortwave radiation.

Fig. 7.

Maps of March–May average of effective white-sky (diffuse) land albedo, model (1999–2008) minus MODIS (2001–06), for (top) LM3.0 and (bottom) LM3.1. Effective albedo is ratio of time averages of up- and downward shortwave radiation.

LM3.0 shows an extensive positive low-latitude bias due to excessive input leaf reflectances. The negative bias in the Sahel stems from excessive vegetation density in LM3.0 at the end of the dry season. Warm-region albedo errors are small in LM3.1, whose leaf reflectances were set on the basis of observed surface reflectances.

d. Interception loss

The differences in interception capacities between LM3.0 and LM3.1 produced corresponding differences in fraction of precipitation lost to interception and subsequent evaporation. In LM3.0, this fraction was 10% for tropical forests, 6% for temperate deciduous forests, and 18% for cold evergreen forests. Except for topical forests, these are considerably lower than observation-based estimates of 10%–15% for tropical forests, 15%–25% for temperate forests, and 25%–35% for cold evergreen forests (Shuttleworth 1988; Dykes 1997; Ward 1975). LM3.1 overall performed better than LM3.0, with average fractions of about 18% in tropical forests, 15% in temperate forests, and 24% in cold evergreen forests.

e. Runoff

River discharge is a natural spatial integral of runoff. As the long-term difference between precipitation and evapotranspiration, it is also an indirect measure of the latter, and, consequently, of the partitioning of net radiation into latent and sensible heat fluxes.

Typical errors in modeled annual-mean runoff (as quantified, e.g., by the basin-area-weighted root-mean-square errors in runoff ratio) in both LM3.0 and LM3.1 are less than 10% of precipitation (Fig. 8). Errors of this size generally are not significant when input precipitation errors, basin discretization, and temporal variability of runoff are taken into account. In particular, with the higher snowfall of LM3.1-S, the overall negative runoff bias (~10% of precipitation) in LM3.1 vanishes; indeed, the 50% increase in snowfall causes a positive bias in the three large cold basins considered. The ~20% negative bias in LM3.1 runoff ratio in one large subtropical basin (that of the Niger River) is explained by a difference in time periods of modeled and observed discharge, together with the large interdecadal change in discharge experienced in the Sahel region; the discharge observations are from the wetter mid-1900s, whereas the model output is for the drier late 1900s.

Fig. 8.

Scatterplots of modeled (1979–2008) against observed (period of stream gauge record) runoff ratio (discharge as a fraction of precipitation for 53 river basins) for LM3.0, LM3.1, and LM3.1-S. Large symbols indicate basins larger than 300 000 km2 in area. Blue symbols indicate cold basins (annual-mean air temperature below 9°C). Red symbols indicate (typically subtropical) basins that are not cold and have high annual precipitation and a pronounced dry season, such that > 0 as defined by Milly and Shmakin (2002). Data are from Milly and Dunne (2002); basin selection began with the set of 82 basins used by Milly and Shmakin, with further limitation to the subset of 53 basins for which the estimated input precipitation error would induce a runoff-ratio error of less than 0.1.

Fig. 8.

Scatterplots of modeled (1979–2008) against observed (period of stream gauge record) runoff ratio (discharge as a fraction of precipitation for 53 river basins) for LM3.0, LM3.1, and LM3.1-S. Large symbols indicate basins larger than 300 000 km2 in area. Blue symbols indicate cold basins (annual-mean air temperature below 9°C). Red symbols indicate (typically subtropical) basins that are not cold and have high annual precipitation and a pronounced dry season, such that > 0 as defined by Milly and Shmakin (2002). Data are from Milly and Dunne (2002); basin selection began with the set of 82 basins used by Milly and Shmakin, with further limitation to the subset of 53 basins for which the estimated input precipitation error would induce a runoff-ratio error of less than 0.1.

The tendency of LM3.0 to produce higher cold-region runoff than LM3.1 (without inflated snowfall) is explained by the rejection of snowmelt infiltration in LM3.0: when snow melts over frozen ground in LM3.0, it immediately produces runoff and streamflow, thereby suppressing later warm-season loss of water to evapotranspiration. In contrast, when snow melts in LM3.1, meltwater infiltrates via macropores into the frozen active layer, where it freezes and is immobilized until thawing reaches its depth a month or two later.

f. Intra-annual variation of storage, streamflow, and evapotranspiration

The combination of precipitation and river-discharge observations does not constrain the timing of evapotranspiration, because changing storage also enters the balance equation. As an estimate of changing storage at spatial scales appropriate for a global model, we use estimates of column-integrated land water storage obtained by satellite gravimetry in the Gravity Recovery and Climate Experiment (GRACE).

The major features in the global pattern of annual range of terrestrial water storage are reproduced well by all LM3 experiments (Fig. 9). Magnitudes of annual range of storage in the tropical rain belt are underestimated a bit by LM3.0, but not by LM3.1; a sensitivity experiment revealed that the main factor in this difference in storage range between configurations was the change in parameters of hydraulic geometry. Higher snowfall of LM3.1-S improves the northwestern North American simulation but causes excessive storage range in much of northern Eurasia, which is relatively well simulated by LM3.0 and LM3.1.

Fig. 9.

Annual range of monthly-mean water storage from GRACE (kg m−2). GRACE observations (University of Texas Center for Space Research, level 2, release 5) were spectrally truncated to degree 60, destriped, and run through a 200-km Gaussian smoother (Landerer and Swenson 2012); model output was filtered in the same way, for consistency with the GRACE data.

Fig. 9.

Annual range of monthly-mean water storage from GRACE (kg m−2). GRACE observations (University of Texas Center for Space Research, level 2, release 5) were spectrally truncated to degree 60, destriped, and run through a 200-km Gaussian smoother (Landerer and Swenson 2012); model output was filtered in the same way, for consistency with the GRACE data.

To look in more detail at the seasonal variations of water balance in LM3, we evaluate the model with respect to monthly streamflow and evapotranspiration, in the context of intra-annual changes of total storage, for six major gauged drainage areas—two each in low, middle, and high latitudes (Fig. 10). For the most part, the observed seasonal phase and amplitude of streamflow, evapotranspiration, and storage change are reproduced very well by LM3.0, LM3.1, and LM3.1-S, albeit with notable exceptions. In the high-latitude Yenisei and Mackenzie basins, LM3.0 produces peak snowmelt-season runoff too early, as a result of runoff above frozen ground, as discussed earlier; consistent with this, modeled storage peaks earlier than observations in these basins. In contrast, the seasonality of discharge is realistic in both LM3.1 and LM3.1-S. The total volume of runoff is much too low in LM3.1, and even deficient in LM3.1-S, suggesting the need for high snow bias correction. However, the storage data (as well as Fig. 9) suggest that LM3.1-S is already overcorrecting for snow bias. Resolution of this discrepancy can be found in the evapotranspiration results, which imply modeled cold regions are losing too much water to the atmosphere during northern summers, thereby reducing runoff below realistic values.

Fig. 10.

Modeled and observed (obs) average annual cycles of water-balance fluxes (mm yr−1) and storage (mm) for six regional-scale drainage areas; obs(S) indicates 50% increased snow input. Precipitation observations (1949–2008) are from input dataset. Discharge (variable period of record) was measured by stream gaging. Storage (2004–08) was measured by GRACE. “Observed” evapotranspiration was computed from observations as precipitation minus discharge minus storage change. Model results are all based on the common time period 2004–08. For precipitation, difference between observation and model is due only to differences in time periods. GRACE observations were scaled to adjust for the effects of filtering, following Landerer and Swenson (2012). Gauge locations are Vicksburg (Mississippi), Manacapuru (Amazon), Brazzaville (Congo), Arctic Red River (Mackenzie), and Igarka (Yenisei). The Intermountain West drainage area is the union of 16 small basins draining mountainous areas of western North America.

Fig. 10.

Modeled and observed (obs) average annual cycles of water-balance fluxes (mm yr−1) and storage (mm) for six regional-scale drainage areas; obs(S) indicates 50% increased snow input. Precipitation observations (1949–2008) are from input dataset. Discharge (variable period of record) was measured by stream gaging. Storage (2004–08) was measured by GRACE. “Observed” evapotranspiration was computed from observations as precipitation minus discharge minus storage change. Model results are all based on the common time period 2004–08. For precipitation, difference between observation and model is due only to differences in time periods. GRACE observations were scaled to adjust for the effects of filtering, following Landerer and Swenson (2012). Gauge locations are Vicksburg (Mississippi), Manacapuru (Amazon), Brazzaville (Congo), Arctic Red River (Mackenzie), and Igarka (Yenisei). The Intermountain West drainage area is the union of 16 small basins draining mountainous areas of western North America.

Discharge can be seen to be a small component of the water balance in the Mississippi River basin, where evapotranspiration almost balances precipitation. Here, small relative errors in precipitation and/or evapotranspiration translate to larger relative errors in discharge; LM3.1-S gives the best simulation of discharge, but the small differences among experiments might not be significant. The amplitude of storage is reproduced well by LM3.1 and LM3.1-S, but is too weak in LM3.0. Evapotranspiration is simulated well. In the Intermountain West, the seasonal cycle of discharge is phased a full 3 months ahead of observations in all experiments. This error appears attributable to insufficient snow storage in the model, which does not resolve high mountains. However, error in the seasonality of snow cover would not explain another common bias in all experiments, which is the failure to produce substantial flows during the low-flow season. The early modeled peak in evapotranspiration is not supported by observations; this is consistent with the hypothesis that the model lacks realistic snow cover in the spring.

In the Amazon River basin, LM3.1 comes much closer to the observed storage variation than LM3.0 does. The weakness of storage in LM3.0 translates to a wide seasonal range in discharge, which is almost in phase with precipitation. Although LM3.1 does better in this regard, it has a semiannual discharge signal that is absent from the observations. The Congo River basin is well simulated, especially by LM3.1, which captures both the amplitude and phase of the semiannual course of discharge.

g. Water-table depths and lake levels

The distribution of modeled water-table depth in LM3.1 is bimodal (Fig. 2); the water table is shallow (~1–3 m) in humid regions and deep (~10–50 m and deepening in LM3.1, ~100 m in LM3.1-DRY) in arid regions. Observational analysis (Fan et al. 2013) confirms the presence of a climatic gradient in water-table depth, albeit relatively subtle, but also reveals much topographically controlled finer-scale variation that is not reproduced by the model. The discrepancy could be partly explained by 1) the coarse grid of the model, 2) the scale mismatch between topographically influenced point measurements and grid-scale, hillslope-average model values, and 3) well-location bias, and 4) the absence of pumping in the model. Nevertheless, the failure of LM3.1 to produce deeper water tables anywhere in humid regions appears to be problematic, as does the sharp gradient in depth between humid and arid regions.

Realism of modeled lake levels varies widely (Fig. 11); r2 values in LM3.1 range from 0 to 0.78 (median 0.30) among the 14 large lakes for which observations were available. The Aral Sea example serves as a reminder of the importance of water management, which is not treated in the model. It should be recalled also that a single stage–discharge relation has been used for all exorheic lakes. Furthermore, as balancers of inflow and outflow, lakes can be sensitive to relatively small fluctuations of inflow, and errors in those fluctuations can be magnified by lake levels.

Fig. 11.

Modeled (LM3.1) and observed time series of water levels in four large lakes. Observations from U.S. Department of Agriculture Global Reservoir and Lake website (www.pecad.fas.usda.gov/cropexplorer/global_reservoir/).

Fig. 11.

Modeled (LM3.1) and observed time series of water levels in four large lakes. Observations from U.S. Department of Agriculture Global Reservoir and Lake website (www.pecad.fas.usda.gov/cropexplorer/global_reservoir/).

h. Thermal state of lakes and permafrost

In LM3.0, which has no background eddy diffusivity to enhance vertical mixing in large lakes, ice cover on the Laurentian Great Lakes is excessive; complete ice cover persists for 3 months in the model, while observations show annual peak ice coverage of only about a third of lake area (Fig. 12). In LM3.1, in which vertical mixing in large lakes has been artificially enhanced by an optimized background eddy diffusivity, the extent and timing of ice cover are reasonably consistent with observations, although ice onset and loss occur about a month too late. LM3.1 also captures well the interannual variations in ice cover on the Great Lakes (Fig. 12).

Fig. 12.

(top) Observed (1973–2012; Wang et al. 2012) and modeled (1979–2008) seasonal cycle of total ice cover on Lakes Superior, Michigan–Huron, St. Clair, Erie, and Ontario. (bottom) Observed and modeled time series of annual-mean ice-cover fraction.

Fig. 12.

(top) Observed (1973–2012; Wang et al. 2012) and modeled (1979–2008) seasonal cycle of total ice cover on Lakes Superior, Michigan–Huron, St. Clair, Erie, and Ontario. (bottom) Observed and modeled time series of annual-mean ice-cover fraction.

Observational analysis shows that the boundary between continuous (>50% coverage) and discontinuous permafrost generally lies on the cold side of the −2°C (the prescribed freezing point of soil water in the model) isotherm of annual-mean air temperature (Fig. 13). This indicates that mean soil temperature is higher than that of the atmosphere, as a result of winter insulation of the ground by snowpack. In LM3.0, however, the permafrost regions extend to the warm side of this isotherm, indicating the soil is colder than the atmosphere. This bias is corrected in LM3.1, which has a lower and more realistic snow thermal conductivity than LM3.0. Consistent with observations, the area of modeled permafrost in LM3.1 is smaller than the area where annual-mean air temperature is less than −2°C, and almost nowhere does permafrost coexist with annual-mean air temperature above −2°C. The insulating effect of the additional snowfall in LM3.1-S causes a small reduction in permafrost extent from that of LM3.1.

Fig. 13.

Observed (Brown et al. 2002) and modeled spatial distributions of permafrost. The variable Ap is the total area of permafrost. Blue line is the −2°C isotherm for mean annual near-surface air temperature (MAAT). Yellow indicates presence of permafrost at MAAT below −2°C, green indicates presence at MAAT above −2°C, and red indicates absence at MAAT below −2°C. Observed permafrost presence determined by classification as continuous, implying greater than 50% area underlain by permafrost. Modeled permafrost presence indicated where soil liquid content is zero for all of 1979–2008 at a depth of 5 m.

Fig. 13.

Observed (Brown et al. 2002) and modeled spatial distributions of permafrost. The variable Ap is the total area of permafrost. Blue line is the −2°C isotherm for mean annual near-surface air temperature (MAAT). Yellow indicates presence of permafrost at MAAT below −2°C, green indicates presence at MAAT above −2°C, and red indicates absence at MAAT below −2°C. Observed permafrost presence determined by classification as continuous, implying greater than 50% area underlain by permafrost. Modeled permafrost presence indicated where soil liquid content is zero for all of 1979–2008 at a depth of 5 m.

In all experiments, permafrost active-layer depths typically reach 0.5–2 m, with minimum values in the coldest regions. Depth of permafrost ranges to a maximum over 60 m in LM3.1; this is considerably less than maximum observed depths as great as 200–400 m, but the difference is attributable to the long-term effect of initialization on deep temperatures in the model (Fig. 1).

5. Discussion

Physically realistic representation of water and energy balances of a large (~5 × 104 km2) grid cell by what is essentially a single-column model is challenging. On one hand, some processes scale well spatially and so are well constrained by relatively well-known column physics and parameters. On the other hand, subgrid heterogeneities, nonlinearities, and unrepresented processes either lead to recalcitrant errors or become aliased into model conceptualizations and parameters that are not well constrained a priori; these are potentially the weak links in the model. In this regard, one of the less robust aspects of LM3’s formulation and implementation is the use of macroporosity to enable both downslope subsurface flow and delay of the snowmelt peak in river discharge. Other examples include suppression of evaporation from soil by litter resistance and enhanced eddy diffusivity of large lakes.

Furthermore, some deficiencies of LM3 behavior are not readily ameliorated by reasonable adjustment of globally constant parameters. The water table in humid regions appears consistently too shallow; this behavior might be related to excessive simplicity of LM3’s representation of the vertical structure of subsurface hydraulic properties, including the (globally constant) shallow depth to bedrock. Vegetation density in the model varies in stepwise fashion between biomes, with consequences for biophysical properties such as albedo; better representation of vegetation-affected processes in LM3 might result from allowance for greater plant functional diversity. The small number of plant functional types in the model might also explain the sensitivity of vegetation to initial conditions.

Errors in model forcing should not be overlooked as potential explanations for systematic errors in model outputs. Credible adjustments of precipitation amounts for observational bias can lead to first-order effects on observable terms in the water balance, impeding detection of model error. Spatial resolution of inputs is another concern; unresolved topographic influences on precipitation and energy balance might account for LM3’s phase error in the seasonal hydrograph of mountainous drainage areas, and finite horizontal resolution limits the fidelity of representation of small river basins.

It is acknowledged that the model evaluation presented here is ad hoc and does not facilitate comparison of model performance with that of other models. Land-model “benchmarking” (Kumar et al. 2012; Luo et al. 2012) is a nascent activity that has potential to enable more systematic evaluative model intercomparisons in the future.

6. Summary

LM3 is a model of terrestrial water, energy, and carbon balances. It includes multilayer representations of temperature, liquid water content, and ice content of snowpack and soil–bedrock; a topography-based parameterization of saturated area and groundwater discharge to streams; transport and storage of runoff in a global river network; and lakes, lake ice, and lake-ice snowpacks that exchange mass and energy with both the atmosphere and the rivers. LM3 dynamically tracks water-table depth, lake levels, and vertical and horizontal extent of permafrost, including active-layer thickness. Temperature of water is tracked throughout the model, and sensible heat of water is consistently included in energy balances. Carbon balance of land-cover vegetation and soil, along with vegetation structure, phenology, and function, are represented as in LM3V that has been described elsewhere.

Experiments with two model configurations—LM3.0 and LM3.1—elucidated model performance and sensitivities. Critical model parameters and features distinguishing performance of the two configurations include leaf reflectance, interception capacity, number of cold months defining limits of cold evergreen forests, snow thermal conductivity, macroporosity, river hydraulic geometry, litter resistance to vapor transport, and background diffusivity of lakes.

In LM3.1, the vertical domain of the model reaches 200 m into the earth. This depth allows direct treatment of the full thickness of the unsaturated zone and the water table in arid regions and, with inclusion of the geothermal heat flux, representation of the full thickness of permafrost in all but the coldest regions. Explicit representation of such a large thickness of crust comes at the cost of very long (≫300yr) model equilibration times, although surface fluxes reach near equilibrium relatively quickly. Spinup-related global water-storage rate after 270 years is sufficient to drive a spurious sea level rise on the order of 0.5 mm yr−1, so care should be taken when using the model in sea level studies.

Experiments with differing initial amounts of vegetation biomass revealed that the modeled forest–grassland boundary at equilibrium depends on initial conditions. The global distribution and seasonal phenology of vegetation in LM3 generally agrees with satellite-derived NDVI, but stepwise spatial transitions between plant functional types cause corresponding transitions in biophysical properties and processes.

Water-balance partitioning in LM3 is sensitive to plant xylem resistance under non-water-stressed conditions, and realistic partitioning is obtained with a realistic input value of the resistance, but only under an assumption of substantial cold-season precipitation observational bias in cold regions. Spatial and seasonal variations in water fluxes in LM3 generally agree with observational estimates from streamflow measurements. Satellite gravimetry shows that the model captures well the amplitude and phase of the annual cycle of total water storage. The systematic high runoff bias of the LaD model in semiarid subtropical basins (Milly and Shmakin 2002) is consistently absent in LM3, which has a more complete treatment of soil-water storage.

Spatial distribution of water-table depth in LM3 follows the climatic moisture gradient as expected, but the water-table depth does not display the smaller-scale spatial variability expected from observations in humid regions, where the model has a shallow bias. Modeled lake level variations are of widely varying quality. Simulation of ice cover on large lakes requires a calibrated parameterization of the unresolved mixing processes. Areal extent of permafrost is sensitive to snow thermal conductivity and is consistent with observations in LM3.1.

Acknowledgments

We acknowledge helpful advice from R. E. Dickinson, D. P. Lettenmaier, S. W. Hostetler, and two anonymous reviewers. Jia Wang provided data on Great Lakes ice cover. Support was provided by U.S. Department of Agriculture, Grant 2011-67003-30373; by the National Oceanic and Atmospheric Administration, U. S. Department of Commerce, Award NA08OAR4320752; and by NASA through the Jet Propulsion Laboratory, Subcontract 1430081.

REFERENCES

REFERENCES
Barraclough
,
G.
, Ed.,
1988
, The Times Atlas of the World. 7th ed. John Bartholomew & Son, 228 pp.
Beven
,
K. J.
,
1986
: Runoff production and flood frequency in catchments of order n: An alternative approach. Scale Problems in Hydrology, V. K. Gupta, I. Rodriguez-Iturbe, and E. F. Wood, Eds., D. Reidel, 107–131.
Beven
,
K. J.
, and
P.
Germann
,
1982
:
Macropores and water flow in soils
.
Water Resour. Res.
,
18
,
1311
1325
, doi:.
Bonan
,
G. B.
,
2002
: Ecological Climatology: Concepts and Applications. Cambridge University Press, 678 pp.
Brown
,
J.
,
O. J.
Ferrians
Jr.
,
J. A.
Heginbottom
, and
E. S.
Melnikov
,
2002
: Circum-Arctic map of permafrost and ground-ice conditions, version 2.0. National Snow and Ice Data Center, Boulder, CO, digital media. [Available online at http://nsidc.org/data/ggd318.]
Brutsaert
,
W.
,
1982
: Evaporation into the Atmosphere. D. Reidel, 299 pp.
Brutsaert
,
W.
,
2005
: Hydrology. Cambridge University Press, 605 pp.
Campbell
,
G. S.
,
1974
:
A simple method for determining unsaturated conductivity from moisture retention data
.
Soil Sci.
,
117
,
311
314
, doi:.
Collatz
,
T.
,
J. T.
Ball
,
C.
Grivet
, and
J. A.
Berry
,
1991
:
Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: A model that includes a laminar boundary layer
.
Agric. For. Meteor.
,
54
,
107
136
, doi:.
Collatz
,
T.
,
M.
Ribas-Carbo
, and
J. A.
Berry
,
1992
:
Coupled photosynthesis–stomatal conductance model for leaves of C4 plants
.
Aust. J. Plant Physiol.
,
19
,
519
538
, doi:.
Cosby
,
B. J.
,
G. M.
Hornberger
,
R. B.
Clapp
, and
T. R.
Ginn
,
1984
:
A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils
.
Water Resour. Res.
,
20
,
682
690
, doi:.
Cuffey
,
K. M.
, and
W. S. B.
Paterson
,
2010
: The Physics of Glaciers, 4th ed. Academic, 704 pp.
Deardorff
,
J. W.
,
1978
:
Efficient prediction of ground surface temperature and moisture, with inclusion of a layer of vegetation
.
J. Geophys. Res.
,
83
,
1889
1903
, doi:.
Delworth
,
T. L.
, and Coauthors
,
2006
:
GFDL’s CM2 global coupled climate models. Part I: Formulation and simulation characteristics
.
J. Climate
,
19
,
643
674
, doi:.
Dickinson
,
R. E.
,
A.
Henderson-Sellers
, and
P. J.
Kennedy
,
1993
: Biosphere–Atmosphere Transfer Scheme (BATS) Version 1e as coupled to the NCAR Community Climate Model. NCAR Tech. Note NCAR-TN-387+STR, 88 pp., doi:.
Dunne
,
J. P.
, and Coauthors
,
2012
:
GFDL’s ESM2 global coupled climate–carbon earth system models. Part I: Physical formulation and baseline simulation characteristics
.
J. Climate
,
25
,
6646
6665
, doi:.
Dykes
,
A. P.
,
1997
:
Rainfall interception from a lowland tropical rainforest in Brunei
.
J. Hydrol.
,
200
,
260
279
, doi:.
Edlefsen
,
N. E.
, and
A. B. C.
Anderson
,
1943
:
Thermodynamics of soil moisture
.
Hilgardia
,
15
,
31
298
.
Emmett
,
W. W.
,
1972
: The hydraulic geometry of some Alaskan streams south of the Yukon River. USGS Open-File Rep.
72
-
108
, 102 pp. [Available online at http://137.229.113.30/webpubs/usgs/of/text/of72-0108.pdf.]
Fan
,
Y.
,
H.
Li
, and
G.
Miguez-Macho
,
2013
:
Global patterns of groundwater table depth
.
Science
,
339
,
940
943
, doi:.
FAO
,
2003
: Digital Soil Map of the World and Derived Soil Properties. Food and Agriculture Organization, CD-ROM, version 3.6.
FAO
,
2006
: Guidelines for soil description. 4th ed. Food and Agricultural Organization of the United Nations, 97 pp. [Available online at ftp://ftp.fao.org/agl/agll/docs/guidel_soil_descr.pdf.]
Farquhar
,
G. D.
,
S.
von Caemmerer
, and
J. A.
Berry
,
1980
:
A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species
.
Planta
,
149
,
78
90
, doi:.
Fekete
,
B. M.
,
C. J.
Vörösmarty
, and
R. B.
Lammers
,
2001
:
Scaling gridded river networks for macroscale hydrology: Development, analysis and control of error
.
Water Resour. Res.
,
37
,
1955
1967
, doi:.
Gardner
,
W. R.
,
1960
:
Dynamic aspects of water availability to plants
.
Soil Sci.
,
89
,
63
73
, doi:.
Garratt
,
J. R.
,
1994
: The Atmospheric Boundary Layer. Cambridge University Press, 316 pp.
GFDL Global Atmospheric Model Development Team
,
2004
:
The new GFDL global atmospheric and land model AM2–LM2: evaluation with prescribed SST simulations
.
J. Climate
,
17
,
4641
4673
, doi:.
Gleeson
,
T.
,
L.
Smith
,
N.
Moosdorf
,
J.
Hartmann
,
H. H.
Dürr
,
A. H.
Manning
,
L. P. H.
van Beek
, and
A. M.
Jellinek
,
2011
:
Mapping permeability over the surface of the Earth
.
Geophys. Res. Lett.
,
38
,
L02401
, doi:.
Henderson
,
F. M.
,
1966
: Open Channel Flow. Macmillan, 522 pp.
Henderson-Sellers
,
B.
,
1985
:
New formulation of eddy diffusion thermocline models
.
Appl. Math. Model.
,
9
,
441
446
, doi:.
Hillel
,
D.
,
1980
: Fundamentals of Soil Physics. Academic Press, 413 pp.
Hostetler
,
S. W.
, and
P. J.
Bartlein
,
1990
:
Simulation of lake evaporation with application to modeling lake level variations of Harney-Malheur Lake, Oregon
.
Water Resour. Res.
,
26
,
2603
2612
, doi:.
Hurtt
,
G. C.
, and Coauthors
,
2011
:
Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands
.
Climatic Change
,
109
,
117
161
, doi:.
Ibbitt
,
R. P.
,
1997
:
Evaluation of optimal channel network and river basin heterogeneity concepts using measured flow and channel properties
.
J. Hydrol.
,
196
,
119
138
, doi:.
Koster
,
R. D.
,
M. J.
Suarez
,
A.
Ducharne
,
M.
Stieglitz
, and
P.
Kumar
,
2000
:
A catchment-based approach to modeling land surface processes in a general circulation model: 1. Model structure
.
J. Geophys. Res.
,
105
,
24 809
24 822
, doi:.
Kumar
,
S. V.
,
C. D.
Peters-Lidard
,
J.
Santanello
,
K.
Harrison
,
Y.
Liu
, and
M.
Shaw
,
2012
:
Land surface Verification Toolkit (LVT)—A generalized framework for land surface model evaluation
.
Geosci. Model Dev.
,
5
,
869
886
, doi:.
Landerer
,
F. W.
, and
S. C.
Swenson
,
2012
:
Accuracy of scaled GRACE terrestrial water storage estimates
.
Water Resour. Res.
,
48
,
W04531
, doi:.
Leopold
,
L. B.
, and
T. J.
Maddock
,
1953
: The hydraulic geometry of stream channels and some physiographic implications. USGS Professional Paper 252, 55 pp. [Available online at http://pubs.usgs.gov/pp/0252/report.pdf.]
Leopold
,
L. B.
,
M. G.
Wolman
, and
J. P.
Miller
,
1964
: Fluvial Processes in Geomorphology. W. H. Freeman and Co., 522 pp.
Leuning
,
R.
,
1995
:
A critical appraisal of a combined stomatal–photosynthesis model for C3 plants
.
Plant Cell Environ.
,
18
,
339
355
, doi:.
Lunardini
,
V.
,
1981
: Heat Transfer in Cold Climates. Van Nostrand Rheinhold, 731 pp.
Luo
,
Y. Q.
, and Coauthors
,
2012
:
A framework for benchmarking land models
.
Biogeosciences
,
9
,
3857
3874
, doi:.
Maherali
,
H.
,
W. T.
Pockman
, and
R. B.
Jackson
,
2004
:
Adaptive variation in the vulnerability of woody plants to xylem cavitation
.
Ecology
,
85
,
2184
2199
, doi:.
Manabe
,
S.
,
1969
:
Climate and the ocean circulation: I. The atmospheric circulation and the hydrology of the earth’s surface
.
Mon. Wea. Rev.
,
97
,
739
774
, doi:.
Meinshausen
,
M.
, and Coauthors
,
2011
:
The RCP greenhouse gas concentrations and their extension from 1765 to 2300
.
Climatic Change
,
109
,
213
241
, doi:.
Miller
,
E. E.
,
1980
: Similitude and scaling of soil-water phenomena. Applications of Soil Physics, D. Hillel, Ed., Academic Press, 300–318.
Miller
,
R. D.
,
1980
: Freezing phenomena in soils. Applications of Soil Physics, D. Hillel, Ed., Academic Press, 254–299.
Milly
,
P. C. D.
, and
K. A.
Dunne
,
1994
:
Sensitivity of the global water cycle to the water-holding capacity of land
.
J. Climate
,
7
,
506
526
, doi:.
Milly
,
P. C. D.
, and
K. A.
Dunne
,
2002
:
Macroscale water fluxes: 1. Quantifying errors in estimation of basin mean precipitation
.
Water Resour. Res.
,
38
,
1205
, doi:.
Milly
,
P. C. D.
, and
A. B.
Shmakin
,
2002
:
Global modeling of land water and energy balances. Part I. The Land Dynamics (LaD) model
.
J. Hydrometeor.
,
3
,
283
299
, doi:.
Pollack
,
H. N.
,
S. R.
Hurter
, and
J. R.
Johnson
,
1993
:
Heat flow from the earth’s interior: Analysis of the global data set
.
Rev. Geophys.
,
31
,
267
280
, doi:.
Schaaf
,
C. B.
, and Coauthors
,
2002
:
First operational BRDF, albedo nadir reflectance products from MODIS
.
Remote Sens. Environ.
,
83
,
135
148
, doi:.
Sheffield
,
J.
,
G.
Goteti
, and
E. F.
Wood
,
2006
:
Development of a 50-yr high-resolution global dataset of meteorological forcings for land surface modeling
.
J. Climate
,
19
,
3088
3111
, doi:.
Shevliakova
,
E.
, and Coauthors
,
2009
:
Carbon cycling under 300 years of land use change: Importance of the secondary vegetation sink
.
Global Biogeochem. Cycles
,
23
,
GB2022
, doi:.
Shuttleworth
,
W. J.
,
1988
:
Evaporation from Amazonian rainforest
.
Proc. Roy. Soc. London
,
B233
,
321
346
, doi:.
Stieglitz
,
M.
,
D.
Rind
,
J.
Famiglietti
, and
C.
Rozenzweig
,
1997
:
An efficient approach to modeling the topographic control of surface hydrology for regional and global climate modeling
.
J. Climate
,
10
,
118
137
, doi:.
Stoy
,
P. C.
, and Coauthors
,
2006
:
Separating the effects of climate and vegetation on evapotranspiration along a successional chronosequence in the southeastern US
.
Global Change Biol.
,
12
, 2115–2135. doi:.
Sturm
,
M.
,
J.
Holmgren
,
M.
König
, and
K.
Morris
,
1997
:
The thermal conductivity of seasonal snow
.
J. Glaciol.
,
43
,
26
41
.
Swenson
,
S. C.
, and
P. C. D.
Milly
,
2006
:
Climate model biases in seasonality of continental water storage revealed by satellite gravimetry
.
Water Resour. Res.
,
42
,
W03201
, doi:.
Tóth
,
J.
,
1962
:
A theory of groundwater motion in small drainage basins in central Alberta
.
J. Geophys. Res.
,
67
,
4375
4387
, doi:.
van der Leeden
,
F.
,
F. L.
Troise
, and
D. K.
Todd
,
1990
: The Water Encyclopedia. 2nd ed. Lewis Publishers, 808 pp.
Wang
,
J.
,
X.
Bai
,
H.
Hu
,
A.
Clites
,
M.
Colton
, and
B.
Lofgren
,
2012
:
Temporal and spatial variability of Great Lakes ice cover, 1973–2010
.
J. Climate
,
25
,
1318
1329
, doi:.
Ward
,
R. C.
,
1975
: Principles of Hydrology. 2nd ed. McGraw-Hill, 367 pp.
Wohl
,
E. E.
,
2004
:
Limits of downstream hydraulic geometry
.
Geology
,
32
,
897
900
, doi:.
Wood
,
E. F.
,
D. P.
Lettenmaier
, and
V. G.
Zartarian
,
1992
:
A land-surface hydrology parameterization with subgrid variability for general circulation models
.
J. Geophys. Res.
,
97
,
2717
2728
, doi:.
Xu
,
J.
,
2004
:
Comparison of hydraulic geometry between sand- and gravel-bed rivers in relation to channel pattern discrimination
.
Earth Surf. Processes Landforms
,
29
,
645
657
, doi:.
Zobler
,
L.
,
1986
: A world soil file for global climate modeling. NASA Tech. Memo. 87802, 33 pp.