A near-global, high resolution land surface parameter dataset for the variable infiltration capacity model

This section describes how we used freely-available data to compile Classic driver input files for the VIC model. First, we created parameter files for VIC-5 Classic, then we converted them to NetCDF format for VIC-5 Image. VIC-5 Classic requires three parameter files: a soil parameter file, a vegetation parameter file, and a vegetation library file. An optional elevation band file can be provided to resolve sub-grid variability in elevation, which is important in regions with complex topography. The parameters are arranged as a relational database: each grid cell has a unique identifier, called a grid cell number, in the soil parameter file, that VIC uses to find the corresponding rows of data in the vegetation parameter and elevation band files. The Image driver uses a different setup, with all parameters stored in a single NetCDF file.

Soil parameters

The soil parameter file for VIC-5 Classic is an ASCII text file that includes soil parameters such as hydraulic conductivity and porosity, but also other kinds of static parameters, such as average precipitation and time zone offset from GMT. Each row of the soil parameter file represents one grid cell, and each column represents a different variable. We compiled the soil parameter file using MERIT²⁰ elevation data, soil texture data from the FAO HWSD, pedotransfer tables relating soil texture to other soil properties, and interpolated weather station data (WorldClim²¹). Any remaining parameters were set to suggested values from the VIC model’s documentation². The following sections describe the estimation of each variable in the soil parameter file, summarized in Table 1.

Elevation and land mask

The VICGlobal soil parameter file uses the Multi-Error-Removed Improved-Terrain (MERIT²⁰) digital elevation model (DEM) to define the elevations, latitudes, and longitudes of each land grid cell. The MERIT DEM is an error-corrected and extended version of the SRTM DEM, with 3 arc-second resolution and coverage from 60°S to 85°N and 180°W to 180°E. Specifically, MERIT is a combination of the SRTM, AW3D, and Viewfinder Panoramas’ DEMs, corrected for striping, speckle, absolute bias, and tree height bias. We used bilinear interpolation to aggregate MERIT to 1/16° resolution and derive a 1/16° MERIT-based land mask and DEM (Figure S1).

Soil texture data

Soil texture (percent sand, silt, and clay) and bulk density were obtained from the FAO HWSD, a gridded soil parameter dataset derived from in-situ measurements of the soil column. We used a 0.05° resolution NetCDF dataset converted from the original HWSD Microsoft Access database by Wieder et al.²². We resampled the HWSD soil data from 0.05° to 1/16° resolution using bilinear interpolation with the MATLAB® function griddedInterpolant. While HWSD has near global coverage, there are missing data in some places around the world, notably Greenland and northern Africa. We filled in these missing data using inpainting, a gap-filling method from the field of image processing. We used the MATLAB® function inpaintnans²³, which uses a partial differential equation method to fill in missing data, to fill gaps in the HWSD data over the MERIT land mask. Figure 1 shows the HWSD bulk density data before and after inpainting.

The HWSD data are divided into “topsoil” and “subsoil” parameters. The first 30 cm of the soil column are considered topsoil and the lower 70 cm subsoil. VIC is typically run with three soil layers, so we created a three-layer soil parameter file by breaking up the 30 cm HWSD topsoil layer into two soil layers: one of 10 cm and one of 20 cm, so the final soil parameter file has three layers, with thicknesses of 10 cm, 20 cm and 70 cm, from top to bottom of the soil column. Ten centimeters has been a common choice for the uppermost layer soil depth in VIC modeling applications since its use by Liang et al.²⁴. Soil layer depths are typically used as calibration parameters. VICGlobal values should be considered a starting estimate.

Calculating soil parameter values based on soil textures

Pedotransfer functions (e.g. Cosby et al.²⁵) relate readily available soil properties, such as soil texture, to less easily-observed properties, such as hydraulic conductivity. After resampling the HWSD data from 1/4° to 1/16° resolution, we estimated soil parameters by classifying each grid cell’s USDA soil texture class and assigning physical soil properties based on a lookup table included with the VIC documentation^2,26. The lookup table (Table 2) relates the 12 USDA soil texture classes to bulk density, field capacity, wilting point, porosity, saturated hydraulic conductivity, and slope of the soil water retention curve in Campbell’s equation. We classified soil textures using the USDA soil texture triangle, as implemented by the MATLAB® function soil_classification²⁷. Figure 2 shows the derived USDA soil texture map. We used these along with the lookup table to estimate saturated hydraulic conductivity (K_sat), the exponent in Campbell’s equation for hydraulic conductivity (expt), fractional soil moisture at the critical point (wcr_fract), where the critical point is about 70% of field capacity, fractional soil moisture at the wilting point (wpwp_fract), quartz content, and porosity for each soil layer. The lookup table²⁶ did not include quartz content, so we supplemented it with the soil texture-quartz content lookup table from Peters-Lidard et al.²⁸.

We set the variable infiltration capacity parameter ({b}_{infilt}=0.2), the maximum baseflow fraction threshold ({d}_{s}=0.001), and maximum soil moisture threshold ({w}_{s}=0.9), their suggested values from the VIC documentation. These parameters, along with maximum baseflow velocity (dsmax) and soil depth, are typically calibrated. We set the baseflow curve exponent c = 2, the soil thermal damping depth dp = 4 m, soil density = 2685 kg/m³, surface roughness = 0.001 m, and snow roughness = 0.0005 m, also based on guidance from the VIC documentation. The soil moisture diffusion parameter phi_s is not used in the current version of VIC, so we set it to the no-data value (−999). The final few soil parameters — dsmax, initial soil moisture (initm), and bubbling pressure (bubble)— were calculated using the following equations, based on guidance from the VIC documentation.

Equation (1) estimates dsmax for each grid cell as the product of soil-column average K_sat and land surface slope, which was calculated from the elevation data using the MATLAB® function gradientm²⁹g. Equation (2), where t_l is the thickness of soil layer l, assumes that initial soil moisture is equal to the fractional soil moisture content at the critical point. Equation (3) calculates bubbling pressure as a function of expt, based on linear regression of bubbling pressure vs. expt³⁰. Figures S2–S9 in the Supplementary Information show maps of each soil parameter. We assumed residual soil moisture, the amount of soil moisture that cannot be removed from the soil by drainage or evapotranspiration, was zero.

Elevation bands

VIC uses an elevation band file (also called a snow band file) to account for subgrid heterogeneity in grid cell elevations. The assumption of uniform elevation over an entire grid cell can lead to modeling errors in mountainous regions, where higher topography is associated with cooler temperatures and higher precipitation rates. The elevation band file accounts for subgrid variability in topography by dividing each grid cell into a number of elevation bands, each of which is simulated separately. VIC adjusts temperature, pressure, and precipitation depending on the elevation in each band. We prepared an elevation band file with five elevation bands by comparing the 1/16° DEM used for the soil parameter file with a 30 arc-second DEM. Both DEMs were derived by aggregating MERIT data. For simplicity, we assumed precipitation was evenly distributed among elevation bands within a grid cell. The elevation band file is provided with the caveat that using elevation bands requires more computing power; users may wish to turn elevation bands on or off (via the VIC global parameter file) depending on their needs.

Vegetation parameters

VIC-5 Classic uses a vegetation parameter file to define the fractional cover of different vegetation types within each grid cell and some of their physical properties. Other vegetation parameters are stored in the “vegetation library” file. (VIC-5 Image simply stores all parameters in a single “parameter” file.) The VIC-5 Classic vegetation parameter file consists of information about fractional cover of each land cover type in each grid cell, and their corresponding root zone depths and root fractions within each root zone. The vegetation parameter file can optionally include time-varying LAI, fractional canopy cover, and albedo data, but it is simpler to specify these in the vegetation library (at the cost of not representing some spatial heterogeneity).

We used MODIS land cover data from the 0.05° MODIS MCD12C1 Collection 6 data product³¹ to assign fractional land cover values to each grid cell by calculating the average land cover for MCD12C1 observations over the 2017 calendar year. We chose 2017 because it was the most recent year with data in all the MODIS-based datasets used for this study, and there is very low interannual variability of land cover³² in MCD12C1 Collection 6. Figure 3 shows majority land cover types from the 2017 MCD12C1 observations.

Like all global land cover data products, MCD12C1 makes classification errors. Sulla-Menashe et al.³² reported 67% overall IGBP classification accuracy for 2001 land cover. Classification errors are more common in the “mixed” land covers, such as cropland/natural vegetation mosaic, shrublands, grasslands, and savannas. Fortunately for our purposes, the vegetation parameters for commonly-confused land covers tend to be fairly similar themselves, which reduces the impact of misclassification on land surface modelling results. For example, the LAI of open shrubland is not too different from the LAI of closed shrubland.

We calculated root fraction as a function of land cover class following the method of Zeng³³, who defined the following formula (Eq. 4) for use in parameterizing land surface models:

where Y = cumulative root fraction, d = depth, and a and b are empirical parameters defined by Zeng³³ for each International Geosphere–Biosphere Programme (IGBP) land cover type, based on a rooting depth database compiled from more than 200 field surveys. We used this formula with depths of 0.1 m, 0.7 m, and dr, corresponding to three root zones. The value of dr, the maximum rooting depth for each IGBP land cover type, was taken from Zeng³³. This method assumes that the depth and distribution of roots depends only on the land cover type; we assume that land cover type is the primary control on root characteristics. Table 3 shows root fractions and root zone depths for each IGBP land cover type.

Like previous large-scale VIC vegetation cover datasets, our vegetation parameter file neglects land cover change over time. However, it does have a few other advantages over past vegetation parameter datasets. The land cover classification used in the N2001 and L2013 VIC parameter sets is referred to as “UMD-NLDAS” because it is a modified version of the AVHRR-based University of Maryland (UMD) land cover product³⁴. The UMD-NLDAS classification was modified for the North American Land Data Assimilation project (NLDAS³⁵) to exclude open water, urban, and snow and ice land cover classes (see BV2019). VICGlobal uses 17 IGBP land cover classes, including urban, barren, perennial snow and ice, and inland water bodies, permitting better description of land cover variability than the 11 UMD-NLDAS classification.

Vegetation library file

The vegetation library maps each land cover type to a set of vegetation parameters (Table 4). We adapted the LDAS vegetation library³⁶ for use with the 17 IGBP land cover classes, taking monthly average LAI, fractional canopy cover (fcanopy), and albedo values obtained from recent MODIS data products. We set architectural resistance (r₀) and minimum stomatal resistance (r_min) to values from literature (described below). The rest of the parameters, which are described in the N2001 paper, were left to their original LDAS vegetation library values. This section describes how we estimated LAI, fcanopy, albedo, r₀, and r_min, and how we transferred the remaining parameters from the 11 UMD-NLDAS land cover classes to the 17 IGBP land cover classes.

We used MODIS observations from the year 2017 to calculate monthly average LAI, fcanopy, and albedo for each IGBP land cover type. We calculated LAI and albedo from the MODIS-based Global LAnd Surface Satellite dataset (GLASS^37,38,39) and fcanopy from NDVI observations (MCD13C1⁴⁰) The expression used for fcanopy follows BV2019:

where NDVI_min and NDVI_max are the minimum and maximum values of NDVI observed for that month. Monthly LAI, fcanopy, and albedo values were calculated by averaging over all grid cells of the same land cover type, counting only cells that were at least 90% homogenous, to avoid noise from grid cells with multiple land covers. Excepting perennial snow and ice land cover, the vegetation parameters in the VIC vegetation library should describe snow-free vegetation. Therefore, before calculating LAI, fcanopy, and albedo for each land cover class, we used fractional snow cover data from MOD10CM⁴¹, a global 0.05 degree monthly snow cover dataset, to exclude grid cells with more than 90% snow cover. Additionally, we set albedo to 0.05 for open water, and we set LAI and fcanopy to 0 for open water and perennial snow and ice.

The resistances r_min and r₀ play a role in determining how much plant transpiration occurs. Higher resistance means less transpiration. Stomatal resistance is resistance to the release of water through the plant stomata, and architectural resistance is the aerodynamic resistance between the leaves and the canopy top⁴². Two sets of resistance parameters have been used in past large-scale VIC implementations. N2001 ran VIC over the entire globe using r_min values adapted from Dorman and Sellers’ global database of r_min values⁴³ computed using the Simple Biosphere Model⁴⁴ (SiB). The Nijssen et al.⁴⁵ r₀ values were taken from Ducoudre et al.’s SECHIBA land surface parameterization⁴². The other set of r_min and r₀ parameters are those used in the LDAS vegetation library and in studies such as Livneh et al.¹². This set of r_min values comes from Mao et al.⁴⁶ and Mao and Cherkauer⁴⁷. We used the r_min values from SiB⁴⁴ and the r₀ values from SECHIBA⁴² for VICGlobal as they appeared to be the better documented values.

For the other parameters in the vegetation library file (displacement height, roughness length, etc.), we assigned values using the existing LDAS vegetation library. Since there are 17 IGBP land cover classes, and only 11 UMD-NLDAS land cover classes in the LDAS vegetation library, we re-assigned some IGBP land cover classes to take the parameters of UMD-NLDAS land cover classes. We remapped barren land, permanent wetlands, snow and ice, urban land, and water bodies to take the parameters of “grasslands” from the LDAS vegetation parameter file. While the characteristics of the barren, snow and ice, urban, and water land cover types clearly differ from those of grasslands, their low LAI and fcanopy values, corresponding to sparse vegetation, essentially “turns off” the other vegetation parameters in the VIC model, as pointed out by BV2019. The other remappings were more straightforward. Croplands and croplands/natural vegetation mosaics inherited values from “croplands,” savannas became “wooded grasslands,” and woody savannas became “woodlands.” We were thus able to assign vegetation parameter values to the each of the 17 IGBP land cover classes.

To calculate global average time series of seasonally-varying vegetation parameters would be of limited interest as the seasonal cycle would average out across the equator. Therefore, we calculated average monthly fcanopy, LAI, and albedo for each vegetation type in each hemisphere, and we developed two separate vegetation library files: one for the northern hemisphere and one for the southern hemisphere. Maps of January and July LAI, fcanopy, and albedo are shown in Fig. 4. For illustrative purposes, the parameter values in this figure have been averaged over the 17 IGBP land cover classes using area-based weighting. Figures S14–S19 show maps of the remaining vegetation parameters. Figures S20–S22 show the cycle of LAI, fractional canopy cover, and albedo for each vegetation type, averaged separately over each hemisphere.

Source: Resources - nature.com