1 Rational

We developed the workflow to be scalable for generating two models of Effective Energy and Mass Transfer (EEMT) (see Rasmussen et al. 2005, Rasmussen and Tabor 2007, Chorover et al. 2011, Rasmussen 2012, Rasmussen et al. 2015).

The two models are called EEMT-Traditional (EEMT-Trad), based on the models given in Rasmussen et al. 2007, Chorover et al. 2011, and Rasmussen 2012; and EEMT-Topographic (EEMT-Topo) (given in Rasmussen et al. 2015)

To calculate the models we are using all open-source GNU licensed software on a distributed High Performance Computing (HPC) framework. The past versions of EEMT were developed using ArcGIS which had some variations in the calculation of solar radiation as well as the climate input data (PRISM vs DAYMET).

2 Locate a Digital Elevation Model within the coterminous USA with DAYMET or PRISM data available.

The first step is to create or download a raster digital elevation models (DEM). Rasters generated by OpenTopography.org are used as examples here, but any properly georeferenced DEM with Coordinate Reference System (CRS) can be used. The climatology data come from DAYMET and are gridded at 1km resolution.

After we have downloaded a sample DEM [suggested resolution is 10 meter (m)] from OpenTopography the script will upscale 1km DAYMET temperature data (min & max) using the reference elevation of the DAYMET pixel and the elevation of the local DEM.

Further calculations include variation in temperature based on topographic position (through aspect insolation/shading effects) and wetness based on the topographic wetness index.

Most of the parameter calculations are done in python scripts (via the GRASS Map Calculator r.mapcalc command).

Attention must be paid to the order of parameters in the r.mapcalc equation line - depending on placement different mathematical operators are executed with different levels of priority. An incorrectly ordered equation will produce the wrong output.

Some parameters require us to generate derivative layers such as slope and aspect from the DEM using other command line tools in either SAGA GIS or GRASS.

3 Matching the Coordinate Reference Systems of both the DEMs and Climatology data

A critical component of developing an interpolation of the climate data involves reprojecting the surface model to the climate observations. Because there is only one DEM, and potentially thousands of climate layers we warp the coordinate reference system (CRS) of the DEM to match the DAYMET data.

Our script uses the gdalwarp command to change the DEM CRS to the same projection system as the DAYMET. Any .TIF file with projection data can be converted to the DAYMET projection using our script.

The DAYMET projection is: 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.

For example:

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

E.g. the script resets the projection from EPSG:26911 (NAD83 UTM Zone 11) to the user-defined Lambert Conformal Conic (lcc) projection that DAYMET uses.

The -tr 10 10 switch sets the pixel size to 10m.

It is important to use the -r bilinear option instead of the nearest neighbor (default).

3.1 Patching files together from OpenTopography

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

4 Locally Corrected Temperature

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

#Locally Corrected Temperature
echo "Calculating locally corrected Temperature"
r.mapcalc "tmin_loc = tmin-0.00649*(dem_10m-dem_1km)"
r.mapcalc "tmax_loc = tmax-0.00649*(dem_10m-dem_1km)"

where:

t_min = temperature minimum for a month from DAYMET. Units are °C, resolution is 1km.
- 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 = high resolution DEM, example set at 10m resolution
dem_1km = DAYMET 1km DEM model (could also be PRISM 800m DEM).
tmin_loc = output file Units are °C resolution is defined by DEM (example used 10m).
tmax_loc = same as for tmin_loc

The input and output files should be in units °C with a typical range of values between -10 and 45 °C for most of North America

Alternatively, we can use the DAYMET 1km product directly - as it has already undergone an elevation-lapse rate atmospheric correction for temperature. If we choose this method we need to first smooth the pixelation which is apparent between 1km pixels such that there are no hard edges in our data. To do this we use a bicubic smoothing spline. We also add a 4km buffer (4 pixel edge) to the tile to make sure the edge effect of the interpolation does not adversely impact the values of our reference DEM surface area.

#Interpolating DAYMET Data to match the DEM: https://grass.osgeo.org/grass64/manuals/r.resamp.interp.html
echo "Interpolating DAYMET Data, increasing region size to include a 4km buffer around DEM"
g.region rast=dem_input -p
#Adds a 4km buffer around the bounding region to make sure enough DAYMET pixels are included for a 16-cell bicubic interpolation over the DEM
g.region n=n+4000 e=e+4000 w=w-4000 s=s-4000 -p
echo "Interpolating TMIN to scale of DEM"
r.resamp.interp tmin output=tmin_int meth=bicubic
echo "Interpolating TMAX to scale of DEM"
r.resamp.interp tmax output=tmax_int meth=bicubic
echo "Interpolating PRCP to scale of DEM"
r.resamp.interp prcp output=prcp_int meth=bicubic
echo "Interpolating VP to scale of DEM"
r.resamp.interp vp output=vp_int meth=bicubic

The calculation for local temperature at this scale instead looks like:

#Locally Corrected Temperature
r.mapcalc "tmin_loc = tmin_int"
r.mapcalc "tmax_loc = tmax_int"

5 Local Potential Evapotranspiration for EEMT-Trad

The equation for daily Potential Evapotranspiration (PET) [mm day^-1] is based on Hamon's eq. Because we are evaluating PET across a monthly time scale these units need to be corrected for the number of days per month to estimate monthly PET.

The locally saturated vapor pressure [kPa] (e_s) parameter must also be calculated prior to generating PET; e_s requires a locally adjusted t_max and t_min to be calculated.

Output value of PET are in mm / day^-1 - need to convert units to mm / month^-1 before the final calculation.

# Local Potential Evapotranspiration for EEMT-Trad
# See Rasmussen et al. 2015 Supplemental 1: https://dl.sciencesocieties.org/publications/vzj/abstracts/0/0/vzj2014.07.0102
# Also see http://www.saecanet.com/20100716/saecanet_calculation_page.php?pagenumber=556
# and http://nest.su.se/mnode/Methods/penman.htm
echo "Calculating local Potential Evaporation-Transpiration (PET) for EEMT-Trad using Hamon's Equation"
r.mapcalc "es_tmin_loc = 0.6108*exp((17.27*tmin_loc)/(tmin_loc+237.3))"
r.mapcalc "es_tmax_loc = 0.6108*exp((17.27*tmax_loc)/(tmax_loc+237.3))"
r.mapcalc "e_s = (es_tmax_loc+es_tmin_loc)/2"
r.mapcalc "tmean_loc = (tmax_loc+tmin_loc)/2" 
r.mapcalc "PET_Trad = (29.8*hours_sun*e_s)/(tmean_loc+273.2)"

where:

es_{tmin_loc} and es_{tmax_loc} are the min and max locally adjusted vapor pressure [kPa].
e_s is the mean saturated vapor pressure e_s [kPa].
hours_sun is the average day length during the month, calculated from the rmean.sh - using the average day length.
PET_Trad must be in units of mm day^-1.

Because PET is equal to zero if air temperatures are negative we can use an 'if' statement to modify the calculation when the input value is less than 0 WITHIN r.mapcalc

6 Calculating Solar Radiation

We calculate solar radiation separately from the EEMT calculation using a script called rsun.sh

To calculate the solar radiation of each pixel we first need to generate slope and aspect of a DEM in decimal degrees using the r.slope.aspect in GRASS or Slope, Aspect, Curvature module in SAGA.

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

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.

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

where

-s is the shadowing effect of terrain turned on

asp_in is aspect input

slope_in is slope input

day is the day of the year

step is the time step, 0.05 is equal to ~3 minutes, 0.25 is 15 minutes, and 0.5 is 30 minutes.

insol_time = the hours per day of sunlight [hours day^-1]

glob_rad = global radiation per day [Wh m² day^-1]

http://whatis.techtarget.com/definition/watt-hour-Wh

The watt-hour is not a standard unit in any formal system, but it is commonly used in electrical applications. An energy expenditure of 1Wh represents 3600 joules (3.600 x 10^3 J). To obtain joules when watt-hours are known, multiply by 3.600 x 10^3. To obtain watt-hours when joules are known, multiply by 2.778 x 10-4

We want the units of R_n to be in Megajoules [MJ], but the units of r.sun are in Watt hours (Wh). So we must convert the units from Watt-hours (Wh) to J and then MJ :

#Local Solar Insolation
echo "Calculating local Solar Insolation in Megajoules [MJ]"
r.mapcalc "total_sun_MJ = total_sun*0.0036"

total_sun_MJ = output converted into megajoules per month

total_sun = average global radiation per day [Wh m² day^-1]

1 Wh = 3600 J
Joules to Megajoules: 3600 J / 1,000,000 J = 0.0036

The rmean.sh can calculate the sum, average, median, standard deviation, and variance in solar radiation across the entire month.

The daily average global radiation [MJ day^-1] is what we use to calculate Penman-Monteith PET.

6.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.2 Locally Corrected Temperature that accounts for Solar Insolation

Correcting local temperature variation based on the aspect of exposure is done by estimating the ratio of obvserved solar radiation versus that of a flat surface. Angles that are pointed toward the sun have greater radiation, and thus a higher air temperature

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:

#Calculate Slope and Aspect
echo "Running r.slope.aspect"
r.slope.aspect elevation=dem slope=slope_dec aspect=aspect_dec
#Create flat map
echo "Creating Flat Map"
r.mapcalc "zeros=if(dem>0,0,null())"
echo "Running r.sun on Flat Map"
#Using dem and flat slope and aspect, generate a global insolation model with local shading off
r.sun elevin=dem aspin=zeros slopein=zeros day=$DAY step=$STEPSIZE dist=$INTERVAL glob_rad=flat_total_sun

#Using dem and slope and aspect (decimal degrees), generate a global insolation model with local shading effects on
echo "Running r.sun using dem, aspect, slope"
r.sun -s elevin=dem aspin=aspect_dec slopein=slope_dec day=$DAY step=$STEPSIZE dist=$INTERVAL insol_time=hours_sun glob_rad=total_sun

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:

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

where

S_i = is a ratio of total shortwave radiation on the observed surface versus shortwave radiation of the flat surface:
The tmin_topo is the same for tmin_loc because minimum temperatures are met at night in the absence of solar influences.

6.3 Calculating Albedo

Rasmussen et al. 2015 use the MODIS MCD43A3 product at 500m.

http://modis-atmos.gsfc.nasa.gov/ALBEDO/index.html

S = (1-a)TrSp

a = albedo (-)

Tr = transmissivity, calculated from ...

6.4 Calculating Linke Atmospheric Turbidity

r.sun contains a default parameter for Linke turbidity, which we have not modified here.

Cloud cover % (monthly) 1km : http://www.earthenv.org/cloud

6.5 Net Longwave outgoign radiation

http://edis.ifas.ufl.edu/pdffiles/ae/ae45900.pdf

L_n = σ * ((t_max+273.16)^4 + (t_min +273.16)^4)/2 * (0.34-0.14 *sqrt(e_a) * (1.35 * (R_s / R_so)-0.35)

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
    - Finds the median TWI of the entire surface.
      #Local Water Balance (accounting for topographic redistribution of run-off) echo "Calculating median of TWI for local water balance" r.mapcalc "twi_median=median(twi)" echo "Calculating the median conservative wetness index" r.mapcalc "a_i = twi/twi_median"
      For further investigation: https://www.researchgate.net/profile/Jonathan_Phillips4/publication/229619411_A_saturation_based_model_for_wetland_identification/links/00b4951e07814409b0000000.pdf

8 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). A nice explanation for calculating Penman-Monteith is here.

The FAO-56 Penman-Monteith method to estimate ET_o can be derived [Eq. 1]:

Where,

ET_o = Evapotranspiration [mm day^-1];

R_n = net radiation at the crop surface [MJ m^-2 d^-1];

G = soil heat flux density [MJ m^-2 d^-1];

T = mean daily air temperature at 2 m height [°C];

u₂ = wind speed at 2 m height [m s^-1];

e_s = saturation vapor pressure [kPa];

e_a = actual vapor pressure [kPa];

δe = es - ea = saturation vapor pressure deficit [kPa];

Δ = slope of the vapor pressure curve, [kPa ºC^-1];

γ= psychrometric constant, [kPa °C^-1];

L = latent heat of vaporization, 2.26 [MJ kg⁻¹];

M = molecular weight of water vapor in dry air, 0.622

P_a = air density [kg m^-3];

C_p = specific heat of moist air [MJ kg^-1 *C^-1];

R_a = s m^-1

The reference evapotranspiration, ETo , provides a standard to which:

• evapotranspiration at different periods of the year or in other regions can be compared;

• evapotranspiration of other crops can be related.

8.1 Raster Calculations

Penman-Monteith incorporates Net radiation (R_n) and vapor pressure as:

PET=(D(R-G)+p_a*c_p*(vp_stopo-vp_loc/r_a))/(L*(m_vp+g_psy)) IMPORTANT: check the units of each of the input terms - we want to try and keep everything in Mega Joules.

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 [J]
G = ground heat flux = 0
p_a = mean air density [UNITS?]
- 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 = temperature in degrees Celsius [C]
  - K = temperature in Kelvin, C +273.125
    echo "Calculating the mean air density" 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 = 0.001013 [MJ kg^-1 C^-1] = 1,013 J kg^-1 C^-1

vp_s_topo = saturated vapor pressure of the topographically modified temperature (different than in section 4 above)

r.echo "Calculating the average temperature corrected vapor saturation"
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*tmin_topo)/(tmin_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
    echo "Calculating the local vapor pressure" r.mapcalc "vp_loc = 6.11*10^(7.5*tmin_topo)/(237.3+tmin_topo)"
r_a = aerodynamic resistance
- r_a = (4.72 * (ln(z_m / z₀))²) / (1+0.536 * U_z)
  - z_m = height of meteorological measurement = 2 m
  - z₀ = 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
      echo "Calculating the aerodynamic resistance" r.mapcalc "ra = (4.72*(log(2/0.00137))*2)/(1+0.536*5)"
L = latent heat of water evaporation = 2.45 [MJ kg^-1] = 2,450,000 J kg^-1
m_vp = slope of the saturated vapor pressure-temperature relationship calculated using mean temperature
- m_vp = 0.04145*exp(0.06088*t_mean)
  - exp is an exponential function for the mathematical constant e = 2.7182818284.
    echo "Calculating the slope of the saturated vapor pressure-temperature relationship" 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 = 0.0065
        local lapse rate, assumed to be 6.5 degrees C km-1
      - z = elevation [m]
  - e = ratio of molar mass of water to dry air = 0.622
    echo "Calculating atmospheric pressure" r.mapcalc "P = 101.3*(293-(0.0065*dem_input)/293)^5.26" echo "Calculating psychromatic constant g_psy [kPa C^-1]" r.mapcalc "g_psy = (0.001013*P)/1.5239"

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

echo "Calculating Actual ET using a Zhang-Budyko curve"
r.mapcalc "AET_zb = if(tmean_topo > 0, prcp*(1+(PET_topo/prcp)-(1+(PET_topo/prcp)^2.63)^(1/2.63)), 0)"

Example code

#Potential Evapotranspiration for EEMT-Topo using Penman-Monteith 
echo "Calculating PET [mm day^-1] for EEMT-Topo using the Penman-Monteith Equation"
echo "Calculating atmospheric pressure"
r.mapcalc "P = 101.3*(293-(0.0065*dem_input)/293)^5.26"
echo "Calculating psychromatic constant g_psy [kPa C^-1]"
r.mapcalc "g_psy = (0.001013*P)/1.5239"
echo "Calculating slope of the saturated vapor pressure-temperature relationship"
r.mapcalc "m_vp = 0.04145*exp(0.06088*tmean_topo)"
echo "Calculating aerodynamic resistance"
r.mapcalc "ra = (4.72*(log(2/0.00137))*2)/(1+0.536*5)"
echo "Calculating local vapor pressure [kPA]"
r.mapcalc "vp_loc = 6.11*(10^(7.5*tmin_topo)/(237.3+tmin_topo))"
echo "Calculating average temperature corrected vapor saturation"
r.mapcalc "f_tmin_topo = if(tmin_topo > 0, 6.108*exp((17.27*tmin_topo)/(tmin_topo+237.3)), 0)"
r.mapcalc "f_tmax_topo = if(tmax_topo > 0, 6.108*exp((17.27*tmin_topo)/(tmin_topo+237.3)), 0)"
r.mapcalc "vp_s_topo = if(tmean_topo > 0, (f_tmax_topo+f_tmin_topo)/2, 0)"
echo "Calculating mean air density"
r.mapcalc "p_a = 101325*exp(-9.80665*0.289644*dem_input/(8.31447*288.15))/(287.35*tmean_topo*273.125)"
echo "Calculating solar influenced PET"
echo "Calculating net radiation in megajoules"
r.mapcalc "total_sun_mj = total_sun*0.036"
r.mapcalc "PET_topo = if(tmean_topo > 0, (m_vp*total_sun_mj+p_a*1013*((vp_s_topo-vp_loc)/ra))/(2450000*(m_vp+g_psy)), 0)"

9 The Final Outputs!

9.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 = (prcp_monthly - PET)*c_w in units of mm month^-1
prcp_monthly = monthly precipitation [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)

#EEMT-Traditional
echo "Calculating the Effective Precipitation term E_ppt for EEMT-Trad"
r.mapcalc "E_ppt_trad = if(prcp > 0, (prcp - (PET_Trad*10))*4185.5, 0)" #PET_Trad was in units of cm/month - converted to mm
echo "Calculating E_bio term for Net Primary Productivity for EEMT-Trad"
r.mapcalc "NPP_trad = if(tmean_loc > 0,3000*(1+exp(1.315-0.119*(tmax_loc+tmin_loc)/2)^-1), 0)"
r.mapcalc "E_bio_trad = if(tmean_loc > 0, NPP_trad*(22*10^6), 0)"
echo "Calculating EEMT-Trad = E_ppt + E_bio"
r.mapcalc "EEMT_Trad = (E_ppt_trad + E_bio_trad)/1000000"
#EEMT-Topographical
echo "Calculating parameters for EEMT-Topo"
echo "Calculating northness"
r.mapcalc "N = cos(aspect_rad)*sin(slope_rad)"
echo "Calculating Net Primary Productivity modified by topography, NPP_topo, based on Whittaker and Niering 1975"
r.mapcalc "NPP_topo = 0.39*dem_10m+346*N-187"
echo "Calculating energy from NPP, E_bio_topo"
r.mapcalc "E_bio_topo = NPP_topo*(22*10^6)"
echo "Calculating effective energy from precipitation, E_ppt_topo"
# At this point everything is still in Joules [J] - convert back to Mega Joules [MJ] 
r.mapcalc "E_ppt_topo = if(prcp > 0, F*4185.5*tmean_topo*E_bio_topo, 0)"
echo "Calculating EEMT_topo"
r.mapcalc "EEMT_Topo = (E_ppt_topo + E_bio_topo)/1000000"

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⁶

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)

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.

9.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
    r.mapcalc F=a_i*prcp
DT = T_ambient - T_ref [K]
- T_ambient = ambient temperature at time of water flux [K]
- T_ref = 273.15 K
```
r.mapcalc DT=((tmax_topo+tmin_topo)/2) - 273.15
```
  If tmax_topo and tmin_topo are in units of C you can remove the T_ref value.
```
 r.mapcalc DT=((tmax_topo+tmin_topo)/2)
```
  Reordering the terms and running them in order:
```
r.mapcalc F=a_i*prcp
r.mapcalc DT=((tmax_topo+tmin_topo)/2)
r.mapcalc E_ppt=F*4185.5*DT*E_bio
```
  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 multiplied by the cosine of aspect (in radians):

N=sin(slope_radians) * cos(aspect_radians).

r.mapcalc N=cos(aspect_rad)*sin(slope_rad)
r.mapcalc NPP_topo=0.39*dem_10m+346*N-187

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

Rasmussen et al. (2015) Quantifying Topographic and Vegetation Effects on the Transfer of Energy and Mass to the Critical Zone. doi:10.2136/vzj2014.07.0102