Creating Soil Parameter File

Note

This document needs to be updated to reflect current practices and data sets.

This page describes methods for developing your own soil parameter file, and provides links to existing soil parameter files.  It is assumed that the reader has completed defining their study domain as described in VIC model Setup (Step 1) - Establishing the Domain before completing the steps in this tutorial.

Steps here assume that the user is working with ArcGIS and has the ArcHydro toolset installed.  Most of the steps involved should translate directly to geoprocessing algorithms available in QGIS, SAGA and other commonly available, non-proprietary GIS packages.

Sources of soil data

This is not a comprehensive list of soil databases, please contribute additional sources by posting tot he page comments.

• Globally there is the Harmonized World Soil Database (HWSD), http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/HWSD_Data.html?sb=4. This replaces the earlier FAO sponsored databases.  Data is gridded to 30 m resolution globally, but that does not mean that the source soils data is at any higher resolution than in previous versions.

  • Citation: FAO/IIASA/ISRIC/ISSCAS/JRC, 2012. Harmonized World Soil Database (version 1.2). FAO, Rome, Italy and IIASA, Laxenburg, Austria.
  • [Soils data manual](https://www.fao.org/fileadmin/templates/nr/documents/HWSD/HWSD_Documentation.pdf}

• Within the United States there are many products produced by the NRCS or that have been produced based on the NRCS database.

  • SSURGO is based on the county level is the finest scale database available

  • Efforts are underway to grid this dataset, so gSSURGO is coming on-line for many states, we have data for the following states stored on the RCAC shared drive:

    • Illinois
    • Indiana
    • Ohio
    • Wisconsin 

  • STATSGO is a coarse resolution state based product. Updates to SSURGO move into the STATSGO database very slowly, so SSURGO is more up to date.

  • CONUS-US is a Penn State effort to produce a gridded data set for hydrology modelling applications. It was a good product at the time, but I don\'t know that it has been updated since it original creation.

  • POLARIS dataset

Creation of a VIC model soil parameter file

This discussion and the provided scripts were developed for the HWSD but should be applicable to any other high-resolution gridded soils database that makes use of common standards.  For instance, I expect that everything provided here would apply to the gSSURGO database with only minimal modification.

The instructions for processing soils data into a VIC model soil parameter file requires four different files:

1. High resolution file with soil IDs (MUID, MU_GLOBAL, etc) for each element
2. High resolution mapping of VIC cell numbers to the same elements used by the soil file
3. A file with cell numbers at VIC resolution
4. A comma separated variable (CSV) file of the soil properties.

The first three files can be provided as ArcInfo ASCII raster grid files, or XYZ files, where X is 'lat', Y is 'lon', and Z is 'value'.  Note that these names must appear in the header line of the input files.  For gridded data the first two files must have the same number of cells, if working with XYZ data then the first two files should have the same number of lines.  The third file should have fewer cells or lines than the first two files.  All three spatial files must be in the same input format, so either all Arcinfo ASCII grids, or all XYZ text files (white space separation, and columns headers of 'log', 'lat' and 'value')

Extraction of soils for the VIC model domain

General instructions for setting up VIC model can be found at VIC model Setup (Step 1) - Establishing the Domain and should be completed before following these steps.

It can be a challenge to get the two high-resolution grids to match up exactly. Here are some strategies that seem to help:

1. Resample/convert the VIC cell number data to a high-resolution raster that is a multiple of the soil dataset.  For example if final simulation resolution is 1 and soil file is 0.45, then select 0.05, or 0.01 (finer resolution will result in more robust grid development, but increase processing time and demand for resources), rather than 0.1.
2. Resample the soil grid to the same high resolution before clipping to the study domain.  Carefully set the extent of processing, and make use of the resampled VIC cell number raster as a snapping grid (both set in the Environment panel under Processing Extent of the tool being used).
3. Once resampling is completed, use the Extract by Mask tool to pull out the high-resolution soil grid for further analysis. 
4. Before proceeding, double check that the extent and the number of rows and columns agree between the high-resolution soil and cell number maps.  If not, you need to repeat these steps using a high resolution.

High resolution mapping unit (MU) grid

The gridded soil database will have a base raster layer containing some form of the Mapping Unit (MU) ID.  The HWSD database uses MU_GLOBAL, while most US databases use the simpler MUID.  This file should be at a higher resolution than the VIC model solution elements, so that multiple MU IDs will appear in each VIC element.  This file should be projected and clipped to match the Analysis Region.  Then it should be exported as an ASCII grid file using the Raster to ASCII conversion tool in ArcGIS. 

If your final mapping units are not gridded, then you should use the arcgrid_to_xyz.scr tool to convert the file into an XYZ file.  Add the headers 'lon', 'lat' and 'value' to the first line of the file.

WARNING: If you need to resample this grid make sure to preserve the MU ID values, do not calculate statistics on the MU ID values as they are indices to the soil database.

High resolution VIC cell number grid

You will need a raster file that exactly matches the extent and resolution of the MU grid, but where the grid values are equal to the VIC cell numbers.

• If using a gridded model domain, use the resample command in ArcGIS to sample the VIC model resolution grid cell number raster to the resolution of the MU ID grid.  If there is not an integer scaling factor between the native resolution of the MU ID grid and the VIC simulation grid, you may need to resample the soil grid to a finer resolution to make sure that the two grids match exactly.
• If using polygons to define the VIC model solution elements, then you can use the Polygon to Raster conversion tool to create a raster grid that matches the extent and cell size of the soil MU ID grid.  WARNING: If using polygons, be sure to use the HUC or other common identifier, not the FID (file index) as some polygon pieces may share a common identifier but be listed as separate polygons.  Using the FID can then result in parts of the HUC or common polygon region being assigned different soils (which VIC cannot handle), or no soil if they are for example islands along a coast line not directly corresponding to a defines soil in the soil database.

Original VIC cell number grid

The resolution of the first two files should be matched so that the Python script can associate MU IDs to grid cells, and determine the mixture of soils within each VIC solution element, without having to do the spatial analysis. This file should have the VIC cell numbers at the actual resolution of the simulation.  This will be used as a template by the Python programs to produce a final output files at the correct resolution.  This file should be either an ASCII raster if the VIC model is being run on a grid, or an XYZ file if the VIC simulation is using a non-gridded domain.  To export polygon centers to an ASCII file, use the "Export Feature Attribute ro ASCII" tool, this produces an output file with XCoord and YCoord as the first two columns (representing polygon centroid location) followed by any attributes selected from within the tool.

Export of the soil database

The final pre-step is to export the soil database to an ASCII table that can be parsed by the Python code.  Soil parameters that are required for processing are:

• The number and depth of defined soil layers

• • For the HWSD there are two layers: Top soil, and subsurface soil.  They are 30 and 70 cm, respectively.   They are consistent across the simulation domain.

• For each defined soil layer, you must provide

• • Percent sand, clay and organic matter.
  • Bulk density
  • [ Check for others in file header ]{style="color: rgb(255,0,0);"}

Export database from ArcGIS and a DBX file.  I exported from ArcGIS to DBX, and opened the DBX file in Access, then used Access to export it to a CSV file.  Double check stages of the format, and what exactly my programs are looking for.  If soil properties are stored in separate database files (gSSURRGO?) then they will need to be merged into a single table where the Mapping Unit (MU ID) is the first column.

Simplification of the soil database

The script SimplifySoilDatabase.py was written to simplify the full database for the harmonized world.  Most soil databases, and especially the more recent gridded versions are capable of representing more than one soil type in a grid cell (even at a 30 m resolution).  In most cases this is a reflection of the probability of a specific soil type being in that grid cell.  This is a reflection of the mapping process, which samples soils based on landscape and land cover features, and then uses those parameters to map the point measurements spatially.  The closer and more similar a grid cell is to an actual measurement location the higher the probability that the cell contains only that soil type.

This one to multiple mapping technique (one MU can contain multiple soil types) poses problems for hydrology models like the VIC model.  To build the VIC model soil parameter file, we need to reduce (or simplify) the soil database so that each VIC model solution unit is represented by a single soil type.  The rule for simplification are as follows:

• Grid MU_GLOBAL values can map to multiple soil types with "SHARE" indicating the fraction of the grid cell represented by the current soil, and "SEQ"  indicating its rank in the grid cell.
• Non-soils, so those with no depth defined are removed.
• Soil types with shallow reporting depths (e.g., less than 100 cm) will be reduced in priority.
• Majority soil type will be assigned to each VIC solution element.

The script CreateSoilSimplifyControlFile.py was written to help in the production of the control file required by the script SimplifySoilDatabase.py.  It should be run in the directory from which SimplifySoilDatabase.py will be run.  It will prompt the user to provide the names of input files, and set parameters specific to the soil database being processed.  As of 18 Mar 2016 only the Harmonized World Soil Database (HWSD) is supported.

Additional parameters required for VIC soil parameter file

The script SimplifySoilDatabase.py can only produce inputs related to soil properties for the VIC model soil parameter file.  Before the next processing step can be completed, you will need to create files matching those output from SimplifySoilDatabase.py that include the following supplementary information:

1. Elevation
   • Used to estimate elevation effects on precipitation and temperature, most important in regions with significant topography where snow bands may be used. 
   • Compute mean elevation for each VIC model simulation cell.
2. Slope
   • This is used for calculating Dsmax, but what should be done as we shift to HUCs, since slope is more of a grid based derivative?  Should this perhaps be the median or mean slope within the VIC simulation cell?  Talk with Bowling, Laura C{.confluence-userlink .user-mention username="bowling" linked-resource-id="25631619" linked-resource-version="1" linked-resource-type="userinfo" base-url="https://wiki.itap.purdue.edu"} about this.
   • Compute slope in percent from a projected (coordinates in distance units not degrees) raster elevation grid.
   • Return to VIC simulation coordinate system and estimate median slope for each VIC model simulation cell.
3. Mean annual soil temperature
   • Optional: Required if running with full energy)
   • Typically estimated from mean annual air temperature, though in regions with significant frozen soil air temperature will underestimate soil temperature.  This is because the latent heat of fusion affects freezing in the soil, especially in wet soils, keeping the soil warmer than the air.  Significant snow cover will also keep soil warm.  So adding a few degrees Celsius to average annual air temperature in cold regions can improve simulation of soil frost.
   • The program [XXXXX]{style="color: rgb(255,0,0);"} can be used to calculate average annual meteorological values required here.
4. Mean annual precipitation
   • Mean annual total precipitation for each VIC model simulation cell.
   • Should be estimated along with mean annual air temperature using [XXXXXX]{style="color: rgb(255,0,0);"}.

5. Hour offset from GMT

   • Offset from Greenwich Mean Time (GMT) in integer hours.

   • Used to adjust sub-daily forcing data to match local solar conditions.

   • This link will take you to an ESRI page to download a global map of time zones{.external-link rel="nofollow"}, which can be used to help develop this layer.

   • When working with only daily forcing data on smaller watersheds (those not cover many time zones), this layer can be set to all '0' as the different from GMT will not affect the simulation.

     Important update in VIC 4.1.2

     VIC ADMIN have adopted the following convention for off_gmt:

     • An off_gmt value of 0 indicates that the model start date/time is relative to Greenwich Mean Time (GMT).
     • An off_gmt value of (grid_cell_longitude*24/360) indicates that the model start date/time is relative to local time.
     • When outputting sub-daily results, VIC's output files are referenced to the model start date/time; therefore they are controlled by off_gmt (off_gmt=0 means VIC results are referenced to GMT; off_gmt=(grid_cell_longitude*24/360) means VIC results are referenced to local time).
     • Daily supplied forcings are assumed to start/end at midnight in local time; the forcing start date/time is thus in local time. When VIC disaggregates daily forcings into sub-daily forcings, off_gmt will be used to determine the time lag between the start of the forcing's diurnal cycle and the start of the VIC simulation.
     • Sub-daily supplied forcings are assumed to occur relative to the time zone indicated by off_gmt. Therefore, if VIC outputs these sub-daily forcings, they will occur at the exact same time of day as in the input files. 

     Therefore, if mixing daily and sub-daily forcing inputs, it is important that any sub-daily forcing inputs be shifted as necessary to be in the time zone indicated by off_gmt.

6. July average air temperature

   1. Optional: Required if estimating treeline within the VIC model, i.e. COMPUTE_TREELINE is set to TRUE.
   2. Calculate from the input forcing data, similar to mean annual precipitation and air temperature layers.

Details for most of these fields can be found at the VIC model documentation page on the soil parameter input file.

All additional parameters MUST match the file format and naming conventions used by the output from SimplifySoilDatabase.py - so they should either be ArcGIS ASCII raster files, and XYZ files for ungridded datasets.

Calculation of soil-water characteristics from basic properties

The script CalculatePerSoilDBLayer.py converts the files extracted from the soil database using the script SimplifySoilDatabase.py into a tale of soil parameters for all activated VIC model simulation cells.  Missing soil parameters will be computed using the functions from Saxton and Rawls (2006) and compiled into the library SoilWaterEquations.py stored in /depot/phig/apps/source/python/lib.  These calculations assume the simpler representation of K_theta (the decrease of water conductivity as soil water becomes less than saturation as described by Campbell (1976).  This method does not make use of residual water, so that value is set to zero for soil parameter files generated using this program.

• The program accepts individual files for each soil parameter, either as ArcInfo ASCII raster files or as XYZ files.  It might be useful to convert it to handle a single input table with all parameters, but then additional data types have to be added as extra columns.  Right now separate files for elevation, annual average temperature and precipitation, July temperature, and GMT offset must be created separately, and identified in the OtherFiles array within the script. 18 Mar 2016 by Keith A. Cherkauer.

  • Organic matter should not exceed 8% by weight as solid with such high contents were not used in developing the equations.
  • As of 18 Mar 2016 I do not know what to do about estimating the bulk density of organic soil.  Will have to look more closely at the VIC model source code (Cherkauer, Keith A{.confluence-userlink .user-mention .current-user-mention username="cherkaue" linked-resource-id="21233665" linked-resource-version="1" linked-resource-type="userinfo" base-url="https://wiki.itap.purdue.edu"})
  • As of 18 Mar 2016 the program must be edited to change the files on which it operates, hopefully this will be changed soon.

The script CreateCalculatePerSoilDBLayerControlFile.py was written to help in the production of the control file required by the script CalculatePerSoilDBLayer.py.  It should be run in the directory from which SimplifySoilDatabase.py was be run, and from which CalculatePerSoilDBLayer.py will be run.  It will prompt the user to provide the names of input files, and set parameters specific to the soil database being processed.  As of 18 Mar 2016 only the Harmonized World Soil Database (HWSD) is supported.

Creation of VIC model soil parameter files

• CreateVicSoilParameterFile.py

  • This program uses output from CalculatePerSoilDBLayer.py, which includes a table of soil parameter values for all VIC model cell coordinates and a table that provides links between variable names and their specific values at different depths in the soil column.  The input file provides soil inforamtion at the native vertical resolution of the soil database.  This script is designed to generate a soil parameter file for the number and position of VIC cell soils.  Most often this program will generate a 3 layer file for a traditional VIC model application, or a 12 layer file for applications of VIC-CropSyst, but the number of soil layers is user defined.

  • Modifications:

    • 20150707 Modified to convert run_cell and gridcel columns to type INT.  This eliminates problems reading the file into VIC. Also added step to drop rows/cells that are not going to be run.  This cleans up the final file, if a mask has been provided for the run_cell variable.      KAC

  • As of 18 Mar 2016 This program runs with the default settings from CalculatePerSoilDBLayer.py, but has to be edited to run with any other setup.

• VicPreprocessingLibrary.py

  • This script file contains libraries and other definitions used by the VIC model preprocessing scripts.

  • Modifications:

    • 20150707 Changed name of GetSoilVariable function to GetVariable to reflect that it is being used by more than just the functions devoted to processing soil data.    KAC

Creation of VIC model vegetation parameter files

This content should probably move to Creating Vegetation Files

Vegetation setup files:

• CreateVegparamFileWithoutPonding.py

  • This script uses ASCII reaster grides of fractional land use coverage in IGBP classification to create a standard VIC model VEG_PARAM file.  Ths VEG_PARAM file does not include wetlands / ponded area, but will include urban area as short grass.

• CreateVegparamFileWithPonding.py

  • This script uses ASCII raster grids of fractional land use coverage in IGBP classification to create a standard VIC model VEG_PARAM file.  Ths VEG_PARAM file will include wetlands / ponded area, but will include urban area as short grass.

  • Things to change:

    • need to read in the ponding fraction file
    • sum up the wetland and open water fractions and check that they are not bigger than the ponding fraction file
    • identify largest vegetation fraction, and assign it to be the wetland fraction (so must appear as the first item in the cell vegetation list)
    • extra area can be assigned to the same veg type later in the file

Other programs being documented:

• ClimateStats.py

  • This script is designed to compute annual statistics for climate forcing data.  It will produce files for each VIC simulation cell that include a time series of the annual metric of interest.  These are passed to a program to calculate the Mann-Kendall tau and slope and complete a statistical test to see if there is a statistically significant trend in the annual metric.

    • INPUT: VIC model ASCII forcing file creted for the Tisza River project
    • OUTPUT: Separate files for each metric calculated

  • WARNING: I messed up the water year, which should start October 1, not September 1, so will have to look back at this in the future.

will convert the files extracted from the HWSD and reduced to the majority soil unit per VIC model grid cell using the script SimplifySoilDatabase.py into a VIC model Soil Parameter file.  Currently this works only for gridded soil data!  Missing soil parameters will be computed using the functions from Saxton and Rawls (2006) and compiled into the library SoilWaterEquations.py stored in /depot/phig/apps/source/python/lib.

Creation of the soil parameter file

Checking the soil parameter file

Existing Soil Parameter Datasets

• Copy information on existing soil datasets from PBWorks 17 Mar 2016