Alan F. Hamlet Dennis P. Lettenmaier JISAO/CSES Climate Impacts Group

1. Macro-Scale Hydrologic Modeling: Conceptual Overview and Introduction to the Variable Infiltration Capacity (VIC) Modeling Software Alan F. Hamlet Dennis P. Lettenmaier JISAO/CSES Climate Impacts Group Dept. of Civil and Environmental Engineering University of Washington

2. Outline of the Talk: • Historic Context • Description of the VIC Modeling Package • Model Evaluation • Practical Considerations ftp://ftp.hydro.washington.edu/pub/hamleaf/IDEAM_workshop

3. Historic Context and Conceptual Overview of the Macro-Scale Approach

4. Macro-scale hydrologic models have their origins in the need to simulate the moisture and energy fluxes at the land surface in Global Climate Models (GCMs) as an aggregated quantity over large spatial scales. Early GCMs did not simulate ocean dynamics in an integrated manner and were of very coarse spatial resolution (typically ~500km x 500km grid cells). Computational constraints were limiting. Simple “bucket” land surface models were typical in early implementations in the late 1960s. Atmosphere Moisture Fluxes: Energy Fluxes: Precipitation Evapotranspiration Runoff Radiative Heat Transfer Sensible and Latent Heat Budget Land Surface

5. As the sophistication of GCMs in simulating the global hydrologic cycle and ocean dynamics has developed through time, the need for more sophisticated land surface schemes has been recognized, and computational constraints have also become less important. The development of more sophisticated macro-scale hydrologic models such as VIC (UW), NOAH (NCEP), and CCM3 (NCAR) have resulted. Atmosphere Ocean Moisture Fluxes: Energy Fluxes: Precipitation Evapotranspiration Runoff Radiative Heat Transfer Sensible and Latent Heat Budget Land Surface

6. Conceptual Approach As model grid cells become larger, transfer of water between cells is dominated by water flowing in river channels, and the spatial variability of some physical drivers becomes less important. In particular, the spatial variability of infiltration, runoff, and baseflow within the cell control volume become less important as spatial scale increases.

7. Differences between macro-scale land surface hydrology models and traditional hydrology models

8. Overview of the VIC Model

9. Development of the VIC Model Liang et al. 1994 2-layer soil vegetation model designed to be dynamically coupled to GCMs or weather models (e.g. at 5 degree lat lon resolution) Parameterized infiltration and base flow schemes Single layer energy balance snow model Physically-based vegetation model including canopy effects Physically-based evaporation based on the Penman/Monteith approach

10. Historic Use of the Model Despite the original conception of the model, until very recently the vast majority of the hydrologic research using the model has implemented the model in an “off-line” configuration. That is, driving data is produced (either from observations or simulations) and the model is run as a stand alone tool often as a “black box” used to interpret the hydrologic implications of the variations in the driving data. Most of the improvements in the model have come about because of the discovery of shortcomings of the model during the course of investigations focused on particular “off line” applications. In the last several years, as computational constraints have been relaxed somewhat and the importance of the land surface state as an important driver of atmospheric circulation and precipitation variability, more attention has been focused on using the tool in a dynamic setting. Precipitation and temperature bias remain difficult elements of fully coupled models to resolve. (I.e. it is often difficult to realize the benefits of an improved land surface scheme if precipitation or temperature in the coupled application are strongly biased for other reasons.)

11. VIC Representation of the Landscape

12. Schematic of VIC Hydrologic Model and Energy Balance Snow Model PNW GB CA CRB 12 km 1/8th Deg. 12 km 1/8th Deg. Snow Model

13. Equal Area Elevation Bands The number of bands is determined by the elevation gradient and a specified interval used in pre-processing (e.g. 1500 m/ 500m in the example). Having determined the number of bands, the bands are forced to have equal area by ranking the pixels in a high resolution DEM and dividing them into groups within the cell boundaries with equal numbers of pixels. Temperature and precipitation are different in each band, but are keyed to the driving data for each cell. In current model implementations the mosaic of vegetation types is identical in each elevation band. High Elevation Band Medium Elevation Band Low Elevation Band

14. Vegetation Characteristics • The model represents a particular vegetation class primarily by: • Canopy albedo • Seasonal Leaf Area Index (LAI)– can be unique for each cell. • Canopy storage (assumed to be a function of LAI) • Characteristic vegetation roughness and displacement height • Stomatal resistance (evaporative resistance associated with transpiration) • Architectural resistance (evaporative resistance related to humidity gradient within the canopy structure as compared to the free air) • Rooting depth • Radiation attenuation factor (used to attenuate incoming solar radiation)

15. Representation of Soil Column True depth and composition of the soil column is usually imperfectly known. Porosity, Ksat, field capacity, wilting threshold, residual capacity and other soil characteristics are determined from estimates of soil composition Storage capacity of each layer is depth times porosity. Rooting distribution is specified in the vegetation file as the fraction of the roots occurring in each depth range. The model then calculates the fraction of roots in each soil layer. Thus the rooting depths and soil layers can be varied independently. ~10cm Infiltration and surface runoff Interflow processes ~20cm Baseflow processes ~1.5 m

16. Model Combinatorial Algorithm Each cell is completely independent of the others. The model solves the water and energy balance independently for each elevation band and vegetation type within the cell (plus bare soil). Band 1 Band 2 . . . Band N Then in each time step the model creates a linear combination of each variable according to the fraction of the cell area that is associated with each band and veg type. Veg 1 . . Veg M Final Model Output Value Veg 1 . . Veg M Area fraction weighting by variable Veg 1 . . Veg M

17. Simulation Modes Water Balance Mode: Assumes the surface temperature is equal to the air temperature and solves the water balance. The snow model, however, is always run as an energy balance computation. Full Energy: Solves the surface energy balance to determine surface temperature. A number of options are available for simulating the subsurface heat budget and ground heat flux algorithms. See: http://www.hydro.washington.edu/Lettenmaier/Models/VIC/Technical_Notes/NOTES_model_modes.html

18. Vegetation Canopy

19. Representation of the Canopy and Canopy Storage Precipitation Canopy evap (wet canopy or snow) Transpiration (dry canopy) Canopy Storage (determined by LAI) Canopy “throughfall” occurs when additional precipitation exceeds the storage capacity of the canopy (rain or snow) in the current time step.

20. Infiltration, Runoff and Base Flow

21. The Variable Infiltration Capacity Curve W1 is determined by the soil depth and porosity. Selecting b determines Imax. W1 = 50 mm b = 0.2 Im = W1 * (1+b) = 60 mm

22. The Variable Infiltration Capacity Curve W1 = 50 mm; B = 0.2 Storm 3 Storm 2 Storm 1

23. The Variable Infiltration Capacity Curve W1 = 50 mm; B = 0.5 Storm 3 Storm 2 Storm 1

24. Three Parameter Non-linear Baseflow Relationship The modeler selects Dmax, Ds, Ws. Wmax is determined by the soil parameters. Ws and Ds determine the x and y positions of the linear threshold of the curve. Dmax determines the maximum base flow when the lower layer is fully saturated. Dmax = 100 Ds = 0.3 Ws = 0.5 Ds * Dmax Ws * 500

25. Three Parameter Non-linear Baseflow Relationship The modeler selects Dmax, Ds, Ws. (Wmax is determined by the soil parameters.) Ws and Ds determine the x and y positions of the linear threshold in the curve. Dmax determines the maximum base flow when the lower layer is fully saturated. Dmax = 100 Ds = 0.2 Ws = 0.8

26. Energy Balance Snow Model http://www.ce.washington.edu/pub/WRS/WRS161.pdf

27. Partitioning of Rain and Snow The model currently uses a very simple partitioning method to determine the initial form of the precipitation. E.g. RainMin= 0.0 C SnowMax = 2.0 C If T <= RainMin then 100% snow. If T >= SnowMax the 100% rain. Values in between are a linear interpolation between the two values. E.g. simulated precipitation at 0.5 degrees C would produce 75% snow, 25% rain.

28. Effects of Forest Canopy on Snow Accumulation Loss of canopy increases the snow water equivalent and increases the rate of melt. Source: Storck, P., 2000, Trees, Snow and Flooding: An Investigation of Forest Canopy Effects on Snow Accumulation and Melt at the Plot and Watershed Scales in the Pacific Northwest, Water Resources Series Technical Report No. 161, Dept of CEE, University of Washington. http://www.ce.washington.edu/pub/WRS/WRS161.pdf

29. Evaluation of the Snow Model for Below Canopy and Shelterwood Areas for a Site in the Cascades Source: Storck, P., 2000, Trees, Snow and Flooding: An Investigation of Forest Canopy Effects on Snow Accumulation and Melt at the Plot and Watershed Scales in the Pacific Northwest, Water Resources Series Technical Report No. 161, Dept of CEE, University of Washington. http://www.ce.washington.edu/pub/WRS/WRS161.pdf

30. Evapotranspiration in VIC model wet canopy evaporation ET dry canopy transpiration bare soil surface evaporation

31. Evaporation and Transpiration Evaporation from wet vegetation and transpiration from dry vegetation are estimated by the physically-based Penman Monteith approach. The equation has the form: Evap = (Term1 + Term 2) / (Term 3) (see e.g. equation 3 in Wigmosta et al. 1994) Term 1 is net radiation term, which is primarily a function of incoming solar radiation (cloudiness) and the slope of the saturated vapor pressure-temperature curve. Term 2 is the vapor pressure deficit term which is primarily a function of the humidity and temperature of the air, scaled by an aerodynamic resistance term related primarily to wind speed and surface roughness. Term 3 is a function of the slope of the saturated vapor pressure and resistance terms associated with canopy resistance and aerodynamic resistance Bare soil calculations are similar but include a resistance term related to the soil’s ability to deliver moisture to the surface (a function of upper layer moisture content and soil characteristics)

32. Overall Modeling Structure for Evaporation Calculations Key drivers such as net radiation budget and wind speed are calculated explicitly for each component of the land surface (canopy, understory, bare soil, and snow surface). Wet or dry vegetation is incorporated by selecting the canopy resistance term (same equation). Overstory Understory Wet Vegetation Dry Vegetation Snow No Snow

33. Model Forcing Data • Sub-daily air temperature (°C) • Surface albedo (fraction) • Atmospheric density (kg/m3) • Precipitation (mm) • Atmospheric pressure (kPa) • Shortwave radiation (W/m2) • Daily maximum temperature (°C) • Daily minimum temperature (°C) • Atmospheric vapor pressure (kPa) • Wind speed (m/s) Below is an example of a 4 column daily forcing file: Pcp Tmax Tmin Wind 6.000 22.560 6.440 3.320 1.775 20.800 4.480 1.260 0.000 25.870 4.360 0.970 0.000 28.470 4.610 1.400 0.000 26.130 8.680 0.880 0.500 25.280 6.860 1.770 ...

34. VIC Met Data Preprocessor The driving data for the model can be explicitly given as a time series, or the model will construct a set of complete forcings from a set of limited daily observations (usually daily precip, tmax, tmin, wind speed) following methods developed by Thornton and Running (1997). Hourly temperature data (needed for the hourly snow model simulations) are reconstructed based on empirical relationships to Tmax and Tmin. Cloudiness and solar radiation attenuation and incoming long wave radiation are estimated via the diurnal temperature range. Dew point temperature is related to daily minimum temperature with a long wave radiation correction.

35. Model Evaluation

36. Streamflow Validation Maurer, E.P., A.W. Wood, J.C. Adam, D.P. Lettenmaier, and B. Nijssen, 2002, A long-term hydrologically-based data set of land surface fluxes and states for the conterminous United States, J. Climate. 15, 3237-3251.

37. Evaluation of Streamflow Simulations of the Colorado River at Lee’s Ferry, AZ

38. Comparison of soil moisture simulations and observations for a site in Illinois. Maurer, E.P., A.W. Wood, J.C. Adam, D.P. Lettenmaier, and B. Nijssen, 2002, A long-term hydrologically-based data set of land surface fluxes and states for the conterminous United States, J. Climate. 15, 3237-3251.

39. Trends in April 1 SWE 1950-1997 Mote P.W.,Hamlet A.F., Clark M.P., Lettenmaier D.P., 2005, Declining mountain snowpack in western North America, BAMS, 86 (1): 39-49

40. 1950-1997 relative trends in April 1 SWE vs DJF temperature Obs VIC Obs VIC Obs VIC Obs VIC

41. Practical Considerations

42. Issues Regarding Implementation • The code is freely available to anyone on the web, and some fairly detailed documentation and general support is available via web pages. • Limited technical support is available since the code is maintained by busy grad students and staff researchers whose primary responsibilities lie elsewhere. A large user community shares experience and solutions to common problems, however, and the approach has been effective at resolving most difficulties encountered by users. • Similarly, improvements in the models have frequently come from the user community. • Successful implementation currently requires considerable GIS experience and strong programming skills in a UNIX environment. The ability to handle large data sets using scripts or compiled code is a must. Similar skills are needed to produce driving data sets for the models based on station data or other resources.

43. Computer Issues • VIC runs cell by cell, and can be very efficiently parallelized by dividing the run into separate runs for sub-groups of cells that together cover the entire area of interest. • VIC is typically run in our group on Pentium or AMD architecture using the UNIX operating system (LINEX). LINEX clusters are also being used frequently, but because the runs are executed cell by cell there is not necessarily a great advantage to doing so. • VIC typically uses about 5 meg of RAM when running and RAM usage does not increase with basin size! Considerable disk storage is required for driving data and output, however, and these are dependent on basin size, output time step, etc. • The GCC C compiler (which is available from GNU for free) is specified in the VIC make file and there is little reason to deviate from this choice. Use of another compiler may work, but requires testing. • The model can be successfully run on MS Windows machines. The easiest way to do this is to install CIGWIN, which emulates the UNIX environment. Many pre-processing and post-processing scripts produced by our group, however, require the C shell, which is not identical to the shell used by CIGWIN.

44. Some Practical Considerations Related to Driving Data and Calibration • The quality of driving data sets frequently controls the quality of the hydrologic model simulations, and may also determine the appropriate spatial resolution of the model. (I.e. bad driving data + high resolution model = high resolution junk) • If observed streamflow or other hydrologic data are available, calibration and bias correction can be effective at removing systematic bias from the simulations. If streamflow is the only output needed, bias correction may be preferable to calibration given the many uncertainties in the driving data and observed streamflow records. • Topographic controls on precipitation are of crucial importance in simulating mountain watersheds, and observed data at high elevation is frequently very limited. Statistical approaches like PRISM (Daly et al. 1994) or meso-scale climate model simulations may provide some useful techniques for resolving these difficulties. • For some kinds of studies involving land surface feedbacks or where no driving data is available, it may be preferable to embed the hydrologic model in a meso-scale climate model to produce fully dynamic simulations of climatic drivers and hydrologic variability. Such approaches are expensive, however.