This box searches only this space. The box at the upper right searches the entire iPlant wiki.

Versions Compared


  • This line was added.
  • This line was removed.
  • Formatting was changed.
Comment: Migrated to Confluence 5.3

1 Rational

In order to calculate the two models of Effective Energy to Mass Transfer (EEMT) (Chorover et al. 2011, Rasmussen and Tabor 2007, Rasmussen et al. 2005, 2011, in press), called EEMT-Traditional (EEMT-Trad) and EEMT-Topographically modified (EEMT-Topo) (Rasmussen et al. in press), we will first need to locate and create required input raster surface models. The available raster surface models come directly from, e.g. digital elevation models (DEM), and the climatological observations come from DAYMET. The rest of the data will be created as outputs of calculations using one or more of the OpenTopo and DAYMET surface models as input.

After we have downloaded a DEM [set to 10 meter (m) resolution] from OpenTopography we will upscale the DAYMET climate observation data (e.g. temperature min & max, precipitation) from the 1kilometer (km) National scale DEM. 

To do this we must apply different equations to each 10m DEM pixel which use as a reference the larger 1km DEM pixel value. Final surface models will in some cases require the outputs of other intermediate layers that are computationally very expensive. 

First, to 'upscale' the 1km climate models onto the smaller 10m DEM pixels we will calculate the difference between the reference elevation pixel (1km DEM) and the test elevation pixel (the 10m DEM) using predicted lapse rates of temperature and pressure.orthn

Next, we will further calculate variation in temperature based on topographic position (through solar insolation/shading effects) and local wetness (through the redistribution of water across the surface).

The calculations are done on the command line using the r.mapcalc command.

Make sure to pay attention to the order you give input parameters in the r.mapcalc equation line - depending on placement different mathematical operators are executed with different levels of priority. An incorrectly formatted equation will produce the wrong output.

Some models will require us to generate new layers such as slope and aspect from the DEM using other command line tools in GRASS.

2 Setting the Coordinate Reference System

Before running any r.mapcalc commands you need to use the gdalwarp command to change the OpenTopo DEM to the same projection system as the DAYMET.

The DAYMET projection is the Lambert Conformal Conic, in the WGS_84 datum, with 1st standard parallel at 25 degrees N, 2nd standard parallel at 60 degrees N (latitude of origin 42.5 degrees); and a central meridian of -100 degrees W.

LiDAR data from OpenTopo will be in different projections and datums depending on the various projects they originated from. It will be easiest to set the LiDAR DEMs to the DAYMET projection as there are many more inputs from DAYMET than OpenTopo.

For example:

Code Block
gdalwarp -overwrite -s_srs EPSG:26911 -t_srs "+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -tr 10 10 -r bilinear -multi -dstnodata 0 -of GTiff J:\ACIC2014_midterm\test_folder\opentopo\dem\dem_utm.tif J:/ACIC2014_midterm/test_folder/opentopo/dem/dem_llc.tif
  • That script sets the projection from EPSG:26911 (NAD83 UTM Zone 11) to the user defined Lambert conformal conic (lcc) with the DAYMET defined lat and long boundaries. 
  • The -tr 10 10 switch sets the pixel size to 10m - there was a problem with resampling to a smaller resolution with a change in projection. 
  • It is important you set the interpolation to 'bilinear' instead of nearest neighbor (default). 

2.1 Patching files together from OpenTopography

When files are downloaded from OpenTopo they may be broken up into multiple files that are 'slices' of the selected area.

An easy way to put these back together (if you choose - another option is to try and operate with them broken up), is to use the r.patch command. 

3 Locally Corrected Temperature

For a simple local min temperature calculation the script would look something like: 

Code Block
r.mapcalc tmin_loc=tmin-0.00649*(dem_10m-dem_1km)
r.mapcalc tmax_loc=tmax-0.00649*(dem_10m-dem_1km)


  • tmin = temperature minimum from DAYMET in units °C at 1km resolution. 
  • 0.00649 is equal to the environmental lapse rate of a stationary atmosphere set at 6.49°C per 1,000 m elevation above sea level (rather than two adiabatic lapse rates: dry = 9.8°C and moist = 5.0°C).
  • dem_10m = OpenTopography hydrologically-corrected DEM at 10m resolution
  • dem_1km = DAYMET 1km DEM model
  • tmin_loc = output file in units °C upscaled to 10m resolution
  • tmax_loc = output file "                                                   "

4 Local Potential Evapotranspiration for EEMT-Trad

The equation for Potential Evapotranspiration (PET) is based on Hamon's eq. which appears simple at first. However the locally saturated vapor pressure (vp_s) parameter must first be calculated prior to generating PET. It also requires a locally adjusted tmax and tmin to be calculated.

Code Block
r.mapcalc f_tmin_loc=6.108*exp((17.27*tmin_loc)/(tmin_loc+273.3))
r.mapcalc f_tmax_loc=6.108*exp((17.27*tmax_loc)/(tmax_loc+273.3))
r.mapcalc vp_s=(f_tmax_loc+f_tmin_loc)/2
r.mapcalc PET=(2.1*((hours_sun/12)^2)*vp_s/((tmax_loc+tmin_loc)/2)


PET is in units of mm day^-1

hours_sun must be in units of '12 hours' - thus it has been divided by 12:

Code Block
r.mapcalc PET=(2.1*((hours_sun/12)^2)*vp_s/((tmax_loc+tmin_loc)/2)

5 Local Solar Insolation

To calculate the solar radiation we will be using two commands: r.slope.aspect and r.sun that require a DEM.

The slope and aspect models from r.slope.aspect script must be created prior to running r.sun.

Code Block
r.slope.aspect elevation=DEM slope=slope aspect=aspect

The ideal time step for calculating the sun angle is <4 minutes where step=0.05 is equivalent to 3 minutes. Notice, this is 10x more computationally expensive than running the time step at the default option of step=0.5 or once every 30 minutes. 

Code Block
r.sun -s elevin=DEM aspin=aspect slopein=slope day="1" step="0.05" dist="1" insol_time=hours_sun glob_rad=total_sun  


  • -s is the shadowing effect of terrain turned on
  • aspin is aspect input
  • slopein is slope input
  • day is the day of the year
  • step is the time step, 0.05 is equal to ~3 minutes
  • insol_time = the hours per day of sunlight [hours day]
  • glob_rad = the total radiation per day [Wh m^2 day]
    • 1 Wh = 3600 J

The units of r.sun are in watts per day. Because we want the units to be in joules we must convert the values from Watt-hours (Wh):

Code Block
r.mapcalc total_sun_joules = total_sun*3600
  • total_sun_joules = output converted into joules per day
  • total_sun = total irradiance calculated by r.sun and given in Wh per day.
  • 3600 = 60 minutes * 60 seconds

5.1 Batch processing r.sun models

Some examples of scripts that calculate monthly averages:

GRASS Wiki with code for Seed Maps and Automation

Monthly Average Solar Maps (R.Sun) in GRASS example scripts

Example of a PYTHON script for GRASS r.sun model

6 Locally Corrected Temperature that accounts for Solar Insolation

Correcting local temperature variation using solar insolation is accomplished by multiplying a function of the total solar radiation incoming to a location versus a flat surface to the locally corrected temperature model of elevation.

Because we must calculate a ratio of solar radiation we are going to run r.sun twice. Once using a zeros raster for both slope and aspect without terrain shadows, and once using the actual slope and aspect with terrain shadowing effects turned on.

To create a zeros raster you can use r.mapcalc and write a conditional statement for the DEM:

Code Block
r.mapcalc zeros=if(dem>0,0,null())

The output will be a raster called zeros where if there is an elevation value greater than zero it gets a zero value, else it is null.

Now, rerun r.sun with the zeros raster set as the slope and aspect:

Code Block
r.sun elevin=DEM aspin=zeros slopein=zeros day="1" step="0.05" dist="1" glob_rad=flat_total_sun  


  • flat_total_sun = the total solar radiation calculated for any pixel with 0 slope and 0 aspect and no shading effect

The final mapcalc calculation will multiply the ratio of total sun to flat surface: 

Code Block
r.mapcalc S_i=total_sun/flat_total_sun
r.mapcalc tmin_topo=tmin_loc+(S_i-(1/S_i))
r.mapcalc tmax_topo=tmax_loc+(S_i-(1/S_i))


  • S_i = is a ratio of total shortwave radiation on the observed surface versus shortwave radiation of the flat surface:

7 Local water balance accounting for topographic water redistribution

A standard approach for calculating topographic controls on water redistribution is the 'Topographic Wetness Index' (Bevin and Kirkby 1979, Tarboton 1991) calculated as:

  • TWI = ln(a/tanB) 
    • TWI = index of pixel wetness 
    • a = catchment area in meters square, calculated using the d-infinity multiple flow direction algorithm (Tarboton 1991) 
    • tanB = the tangent of the slope (in degree units)

We can use the TWI model from OpenTopography (produced by the optimized TauDEM) to calculate topographic water redistribution, or we can generate our own model in GRASS with the ta_hydrology module installed in QGIS.

Remember, if you use TWI from OpenTopography you will still need to warp the projection to match the warped DEM and DAYMET data.

The wetting of a pixel is a function of the precipitation 'prcp' minus surface runoff. We will estimate the mass flux, F, of water based on the topographic wetness index, normalized to area.

  • W= AET + F [m s^-1]
  • We can think of total wetting as the actual evapotranspiration (AET) plus F (the subsurface wetting)
    • F = a_i * prcp
      • a_i = mass conservative wetness index (MCWI)
        • a_i = TWI_i / (1/N * sum(TWI_i))  
        • TWI_i = the pixel of interest
        • N is the total number of pixels in the catchment
Code Block
r.mapcalc a_i=TWI/((max(TWI)+min(TWI))/2)

7.1 Potential Evapotranspiration for EEMT-Topo

Rasmussen et al. (in press) calculated the PET rate differently for EEMT-Trad and EEMT-Topo.For EEMT-Topo the PET equation is from Penman-Monteith (Shuttleworth 1993)) incorporates net radiation and vapor pressure as:

  • PET=(D(R-G)+p_a*c_p*(vp_s_topo-vp_loc/ra))/(L*(m_vp+g_psy)) 
Code Block
r.mapcalc g_psy=0.001013*(101.3*((293-0.00649*dem_10m)/293)^5.26)/(0.622*2.45)
r.mapcalc m_vp=0.04145*exp(0.06088*(tmax_topo+tmin_topo/2))
r.mapcalc ra=(4.72*(ln(2/0.00137))2)/(1+0.536*5)
r.mapcalc vp_loc=6.11*10(7.5*tmin_topo)/(237.3+tmin_topo)
r.mapcalc f_tmin_topo=6.108*exp((17.27*tmin_topo)/(tmin_topo+237.3))
r.mapcalc f_tmax_topo=6.108*exp((17.27*tmax_topo)/(tmax_topo+237.3))
r.mapcalc vp_s_topo=(f_tmax_topo+f_tmin_topo)/2 
r.mapcalc p_a=101325*exp(-9.80665*0.289644*dem_10m/(8.31447*288.15))/287.35*((tmax_topo+tmin_topo/2)273.125)
r.mapcalc PET=total_sun_joules+p_a*0.001013*(vp_s_topo-vp_loc)/ra))/(2.45*(m_vp+g_psy)
r.mapcalc AET=prcp*(1+PET/prcp(1(PET/prcp)2.63)(1/2.63))


  • R = topographically modified radiation from r.sun converted into joules
    • from r.sun the units are given in watt-hours  (Wh) per meter square per day (m^-2 d^-1)\
      • 1 Wh = 3600 joules
Code Block
r.mapcalc total_sun_joules = total_sun*3600
  • G = ground heat flux = 0
  • p_a = mean air density
    • p_a = p / R_spec * T
      • p = atmospheric pressure can be derived using different equations 
        • p = 101325 * (1- (0.00649 * elev / 288.15))^(9.80665*0.0289644 / 8.31447*0.00649))
        • alternatively, p ~= 101325 * exp(-9.80665*0.0289644*elev / 8.31447 * 288.15)
      • R_spec = specific gas constant for dry air = 287.35 [J/(kg*K)]
      • T = in degrees K = deg C +273.125
Code Block
r.mapcalc p_a=101325*exp(-9.80665*0.289644*dem_10m/(8.31447*288.15))/287.35*((tmax_topo+tmin_topo/2)+273.125)
  • c_p = specific heat of moist air [MJ kg^-1 C^-1]= 0.001013
  • vp_s_topo = saturated vapor pressure of the topographically modified temperature (different than in section 4 above)
Code Block
r.mapcalc f_tmin_topo=6.108*exp((17.27*tmin_topo)/(tmin_topo+237.3))
r.mapcalc f_tmax_topo=6.108*exp((17.27*tmax_topo)/(tmax_topo+237.3))
r.mapcalc vp_s_topo=(f_tmax_topo+f_tmin_topo)/2 
  • p_loc = local area vapor pressure
    • vp_loc = 6.11 * 10^((7.5*T_d) / (237.3+T_d))
      • T_d = dew point temperature
      • the simple assumption in non-arid sites is: T_d = tmin
Code Block
r.mapcalc vp_loc=6.11*10^(7.5*tmin_topo)/(237.3+tmin_topo)
  • ra = aerodynamic resistance
    • ra = (4.72 * (ln(z_m / z_0))^2) / (1+0.536 * U_z)
      • z_m = height of meteorological measurement = 2 m
      • z_0 = aerodynamic roughness of an open water surface = 0.00137 [m]
      • U_z = windspeed
        • Because we do not have wind speed observations we will assume a constant of 5 m/s
Code Block
r.mapcalc ra=(4.72*(ln(2/0.00137))^2)/(1+0.536*5)
  • L = latent heat of water evaporation = 2.45 [MJ kg^-1]
  • m_vp = slope of the saturated vapor pressure-temperature relationship calculated using mean temperature
    • m_vp = 0.04145*exp(0.06088*tmean)
      • exp is an exponential function for the mathematical constant e = 2.7182818284.
Code Block
r.mapcalc m_vp=0.04145*exp(0.06088*(tmax_topo+tmin_topo/2))
  • g_psy = psychromatic constant 
    • g_psy = c_p*P / e*L
      • c_p = specific heat of moist air (see equation above)
      • P = atmospheric pressure [kPa]
        • P = 101.3 * ((293-n*z) / 293)^5.26
          • n = local lapse rate, assumed to be a constant 6.49 degrees C
          • z = elevation [m] 
      • e = ratio of molar mass of water to dry air = 0.622
Code Block
r.mapcalc g_psy=0.001013*(101.3*((293-0.00649*dem_10m)/293)^5.26)/(0.622*2.45)

The actual evapotranspiration (AET) was described using a Budyko curve (Budyko 1974) that partitions PET and AET relative to an aridity index (ratio of annual PET to annual rainfall). 

  • AET=prcp*(1+PET/prcp-(1+(PET/prcp)^w)^1/w)


  • PET = (see equations above)
  • prcp = precipitation
  • w is an an empirical constant that varies between 2.15 and 3.75, here it is set to w = 2.63 (Zhang-Budyko, Zhang et al .2004)
Code Block
r.mapcalc AET=prcp*(1+PET/prcp-(1+(PET/prcp)^2.63)^(1/2.63))

8 The Final Outputs!

8.1 EEMT-Traditional

EEMT-Trad is reported on in Rasmussen et al. (2011) and defined as:

  • EEMT_Trad = E_ppt + E_bio

where the two terms: E_ppt and E_bio are defined as:

  • E_ppt = (monthly_prcp - PET)*c_w
  • monthly_prcp = total precipitation in mm 
  • c_w = specific heat of water [J kg^-1 K^-1] 
  • 1 calorie/gram °C = 4.186 joule/gram °C 
  • Water (liquid): c_w= 4185.5 J/(kg·K) (15 °C, 101.325 kPa)
Code Block
r.mapcalc f_tmin_loc=6.108*exp((17.27*tmin_loc)/(tmin_loc+273.3))
r.mapcalc f_tmax_loc=6.108*exp((17.27*tmax_loc)/(tmax_loc+273.3))
r.mapcalc vp_s=(f_tmax_loc+f_tmin_loc)/2
r.mapcalc PET=(2.1*((hours_sun/12)^2)*vp_s/((tmax_loc+tmin_loc)/2)
r.mapcalc E_ppt=monthly_prcp - PET

E_bio = NPP * h_bio


  • NPP = mass flux carbon from net primary productivity [kg m^-2 s^-1] = 3000*(1+exp(1.315-0.119*T_mean))^-1
  • h_bio = specific biomass enthalpy [J kg^-1] = 22*10^6 
Code Block
r.mapcalc NPP_trad=3000*(1+exp(1.315-0.119*(tmax_loc+tmin_loc)/2)^-1)
r.mapcalc E_bio=NPP_trad*(22*10^6)


Because EEMT-Trad was only ever calculated on a monthly time-step the iteration down to a daily evaluation of available water from precipitation minus evapotranspiration will need to be somewhat different going forward.

Conceptually, the best way to model this term will be to carry over any positive soil water values from one day to the next after accounting for the previous day's evapotranspiration. This will also eventually require the integration of Snow Water Equivalent (SWE) from negative degree days over the winter leading up to spring snow melt.

8.2 EEMT-Topographic

EEMT-Topo is calculated similarly to EEMT-Trad where: 

  • EEMT_Topo = E_ppt + E_bio

The difference is the way E_ppt and E_bio are derived. For EEMT-Topo the E_ppt term uses a mass conservative wetness index that accounts for the topographic redistribution of water:

E_ppt = F * c_w *DT


  • F = mass flux of water [kg m^-2 s^-1]
    • The F term is derived from the Penman-Budyko approach with the mass conservative wetness index (Section 7.1)
    • F = a_i * prcp
      • a_i = mass conservative wetness index (see equation above in section 6)
      • prcp = daily/monthly precipitation from DAYMET in units of mm
        Code Block
        r.mapcalc F=a_i*prcp
  • DT = T_ambient - T_ref [K]
    • T_ambient = ambient temperature at time of water flux
    • T_ref = 273.15 K
      Code Block
      r.mapcalc DT=((tmax_topo+tmin_topo)/2) - 273.15
      Reordering the terms and running them in order:
      Code Block
      r.mapcalc F=a_i*prcp
      r.mapcalc DT=((tmax_topo+tmin_topo)/2) - 273.15
      r.mapcalc E_ppt=F*4185.5*DT*E_bio
      NPP for EEMT-Topo is calculated using a 'northness' term and elevation:


  • N = Northness is the sine of the slope (in degrees) multiplied by the cosine of aspect (in radians):

N=cos(slope*0.0174532925) * sin(aspect*0.0174532925).

Code Block
r.mapcalc N=cos(slope*0.0174532925)*sin(aspect*0.0174532925)
r.mapcalc NPP_topo=0.39*dem_10m+346*N-187


Chorover, J., Troch, P. A., Rasmussen, C., Brooks, P. D., Pelletier, J. D., Breshears, D. D., ... & Durcik, M. (2011). How water, carbon, and energy drive critical zone evolution: The Jemez--Santa Catalina Critical Zone Observatory.Vadose Zone Journal10(3), 884-899. doi:10.2136/vzj2010.0132

Rasmussen, C., & Tabor, N. J. (2007). Applying a quantitative pedogenic energy model across a range of environmental gradients. Soil Science Society of America Journal71(6), 1719-1729. doi:10.2136/sssaj2007.0051

Rasmussen, C., Troch, P. A., Chorover, J., Brooks, P., Pelletier, J., & Huxman, T. E. (2011). An open system framework for integrating critical zone structure and function. Biogeochemistry102(1-3), 15-29. doi: 10.1007/s10533-010-9476-8