**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 OpenTopography.org, 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:

- 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:

where:

- 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.

where:

PET is in units of mm day^-1

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

**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.

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.

where

- -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):

- 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*

*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:

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:

where

- 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:

where

- 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

- a_i = TWI_i / (1/N * sum(TWI_i))

- a_i = mass conservative wetness index (MCWI)

- F = a_i * prcp

### 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))

where

- 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

- from r.sun the units are given in watt-hours (Wh) per meter square per day (m^-2 d^-1)\

- 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

- p = atmospheric pressure can be derived using different equations

- p_a = p / R_spec * T

- 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)

- 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

- vp_loc = 6.11 * 10^((7.5*T_d) / (237.3+T_d))

- 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

- ra = (4.72 * (ln(z_m / z_0))^2) / (1+0.536 * U_z)

- 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.

- exp is an exponential function for the mathematical constant

- m_vp = 0.04145*exp(0.06088*tmean)

- 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]

- P = 101.3 * ((293-n*z) / 293)^5.26
- e = ratio of molar mass of water to dry air = 0.622

- g_psy = c_p*P / e*L

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)

where

- 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)

**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)

E_bio = NPP * h_bio

where

- 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

#### Comment

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

where

- 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

- DT = T_ambient - T_ref [K]
- T_ambient = ambient temperature at time of water flux
- T_ref = 273.15 K
Reordering the terms and running them in order:NPP for EEMT-Topo is calculated using a 'northness' term and elevation:

NPP=0.39*elev+346*N-187

- 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).

**References:**

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 Journal*, *10*(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 Journal*, *71*(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. *Biogeochemistry*, *102*(1-3), 15-29. doi: 10.1007/s10533-010-9476-8

## 4 Comments

## Daniel Maguire

Has anyone actually been able to get GRASS running on the HPC or India? That seems to be the single biggest holdup right now as we have no way of transitioning our work to something scalable without that piece.

## Tyson Swetnam

Not sure how helpful this will be - but here is an article on using HPC and GRASS: http://link.springer.com/article/10.1007/s12145-014-0165-3

## Weifeng Li

How do we calculate NPP in E_bio?

## Tyson Swetnam

That's awesome you're ready to calculate that!

I've updated the equation with an example of how to calculate NPP.