Surface Renewal Analysis for Energy Flux Exchanges in an Ecosystem: 1: Calculating Ramp Characteristics using R.

July 13, 2017
By

(This article was first published on R-posts.com, and kindly contributed to R-bloggers)

Summary: The application of Surface Renewal (SR) analysis as an alternative to energy and gaseous flux measurements using eddy covariance has gained prominence as a less costly and reliable method.  The R code provided in this post is the first of two posts. In this first post  R programming language is applied to calculate ramp characteristics; namely: ramp amplitude (A) and ramp duration (τ). A sample dataset is provided to demonstrate the functionality of the code.

INTRODUCTION

The surface energy balance equation explains the net solar radiation from the sun to and from the earth’s surface and is comprised of three important elements namely: latent heat (LE), sensible heat (H) and ground heat flux (G). LE is the latent heat flux that results from both evapotranspiration from the soil and vegetation.  Sensible heat warms the layer above the surface while ground flux is associated with heat stored in the soil or pedozone. The sum of these 3 components totals the radiative flux which is a combination of both longwave and shortwave radiation energy. Quantification of these components within a watershed or field, assists scientific understanding of plant development, precipitation and water demands by competing users. There are several micrometeorological methods available to measure these components both directly and indirectly. The eddy covariance (EC) method utilizes an omnidirectional 3D sonic anemometer as well as either an open or closed-path infrared H2O gas analyzing system to measure LE. Solar radiation can be measured using radiometers. Soil heat plates placed at various depths, measure ground heat flux (G) while sensible heat is measured indirectly as the remainder after subtracting G and H from net solar radiation.

The EC methodology is rather expensive and has several limitations which are discussed by several micrometeorologists (e.g. Castellvi et al., (2006); Castellvi et al, (2009a, 2009b); and Suvočarev et al., 2014)). The surface renewal method (SR) however provides an efficient and cheap option for estimating H,  LE energy and methane (CH4) and carbon dioxide (CO2) gas exchange fluxes to and from varied ecosystems such as uniform pine forest (Katul et al., 1996); grapevine canopies (Spano et al., 2000);  rice fields (Castellvi et al., 2006); and rangelands (Castellvi et al., 2008). The method is both cheaper and simple to implement on medium spatial scales such as livestock paddocks and feedlots. SR method does not require instrument levelling and has no limitations in orientation (Suvočarev et al., 2014)).  Additionally, the SR method has demonstrated better energy balance closure (Castellvi, et al., 2008).

To estimate H flux, measurements are taken at one height together with horizontal wind speeds and temperature. To estimate LE fluxes, measurements may be taken at one height together with horizontal wind speeds and moisture concentrations. The Surface Renewal (SR) process is derived from understandings of gaseous flux exchanges to or from an air parcel traveling at a known or specified height over a canopy.  In the SR process an air parcel moves into the lower reaches of the canopy and remains for a duration τ (tau); before another parcel replaces it ejecting it upwards. Paw et al, (1995) created a diagram of this process and likened it to a ramp-like event.

 


Figure 1: Air parcel diagram of surface renewal process.

Figure 2: An ideal depiction of a surface renewal ramp where atmospheric conditions are unstable.

The ramp is characterized by both an amplitude (A) and period (τ – Greek letter :tau). The purpose of this R code is to calculate these two characteristics for measuring both gaseous and energy flux exchanges occurring during the SR process. These values are then applied in estimating H, LE, CO2 and CH4.

Ramp Characteristics

Amplitude

Van Atta (1977) and Antonia and Van Atta (1978) (both cited in Synder et al., 2007) pioneered the surface renewal method. The researchers studied structure functions, ramp patterns and turbulence and their relationships. They utilized the structure function approach to estimate both amplitude and duration from moments recorded in the field (Equation 1). 

(Equation 1)

Where h(V-Vi-j)  is the difference between two sequential high-frequency measurements (time lag interval,  j) of water vapor density or air temperature. M represents the total number of data points in the 0.5 hour time period. i represents the summation index while n is the structure function order. The R code builds of these procedures and calculates the 2nd, 3rd and 5th moments of a structure function following Van Atta (1977).

Applying the three moments to calculate p and q; as:

    (Equation 2)

And;

             (Equation 3)

A cubic equation encompassing ramp amplitude is provided as:

              (Equation 4)

Equating y to 0 and solving for the real roots of the equation yields the ramp amplitude.


Ramp duration

The “inverse ramp frequency” also known as the mean ramp event, τ ,  is calculated as:

(Equation 5)

j represents the instance when the ratio between the third moment and observation time (seconds) is a maximum.


Demonstration

The sample dataset (test.csv) utilized in this package consists of water vapor measurements (gm-3) recorded over 30 minutes (half an hour) at 20 Hz per second. Therefore a total of 36000 measurements are provided. The sample dataset can be downloaded here.


R code

The code provided below converts the H2O concentrations from g per m3 to mmol per m3.

Thereafter, 2nd, 3rd and 5th moments are calculated for each 1/20th of a second. These values are then applied to calculating two values; p and qA (amplitude) is derived by solving for the true roots of a cubic polynomial equation. 

The R code generates a numeric vector that consists of three values:

1) r value – this is the time (seconds) when the ratio between the third moment and observation time (seconds) is a maximum.

2) τ (Greek letter – tau) is the total ramp duration

3) A is the amplitude of the ramp that models the SR process.

In order for the package to function, two supporting packages are necessary. These are the NISTunits and polynom packages.

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, “Package”])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

### usage of packages

packages <- c(“NISTunits”,“polynom”,“bigmemory”)
ipak(packages) ## Loading required package: NISTunits ## Loading required package: polynom ## Loading required package: bigmemory ## Loading required package: bigmemory.sri ## NISTunits   polynom bigmemory 
##      TRUE      TRUE      TRUE require(NISTunits)


 ### cube root function
Math.cbrt <- function(x) 
  {
  sign(x) * abs(x)^(1/3)
  }

####
#### set the directory and read file. In my case my dataset is stored in
#### D:SurfaceRENEWAL folder

setwd(“D:/SurfaceRENEWAL/”)
df<-read.csv(“test.csv”, header=TRUE)

## creating empty data frame
### based on the value of hertz. In this case it is 20 Hz per second..

hertz<-20


dat <- data.frame(col1=numeric(hertz), col2=numeric(hertz), col3=numeric(hertz),stringsAsFactors = FALSE)

####convert to mmol m-3

df$H2O<-1000*df$H2O/18.01528
 
### Creating and Initializing the table with three columns 
dat <- data.frame(col1=numeric(hertz), col2=numeric(hertz), col3=numeric(hertz),stringsAsFactors = FALSE)
dat ##    col1 col2 col3
## 1     0    0    0
## 2     0    0    0
## 3     0    0    0
## 4     0    0    0
## 5     0    0    0
## 6     0    0    0
## 7     0    0    0
## 8     0    0    0
## 9     0    0    0
## 10    0    0    0
## 11    0    0    0
## 12    0    0    0
## 13    0    0    0
## 14    0    0    0
## 15    0    0    0
## 16    0    0    0
## 17    0    0    0
## 18    0    0    0
## 19    0    0    0
## 20    0    0    0 for (p in 1:hertz)
{
  #### initializing cumulative totals
  tempTot2<-0
  tempTot3<-0
  tempTot5<-0
  #### the first column is the frequency
  #### the second column is the duration (in seconds)
  dat$col1[p]<-p
  dat$col2[p]<-p/hertz

  #### initializing temporary variables for calculating 
  #### 2nd, 3rd and 5th moments
  
  tem2<-0
  tem3<-0
  tem5<-0
    for (i in (p+1):nrow(df))
      {
      tem2<-(df$H2O[i]-df$H2O[i-p])^2
      tempTot2<-tem2+tempTot2  
      tem3<-(df$H2O[i]-df$H2O[i-p])^3
      tempTot3<-tem3+tempTot3
      tem5<-(df$H2O[i]-df$H2O[i-p])^5
      tempTot5<-tem5+tempTot5
      }

  dat$col3[p]<-(tempTot2/(nrow(df)-dat$col1[p]))
  dat$col4[p]<-(tempTot3/(nrow(df)-dat$col1[p]))
  dat$col5[p]<-(tempTot5/(nrow(df)-dat$col1[p]))
  
}


dat$moment2<-abs(dat$col3)
dat$moment3<-abs(dat$col4)
dat$moment5<-abs(dat$col5)

#### Solution following Synder et al.(2007)
#### values for p and q
dat$p<-10*dat$col3 – (dat$col5/dat$col4)
dat$q <- 10*dat$col4
#### solutions to the cubic equation after Synder et al., (2007)
dat$D<-(((dat$q)^2)/4)+(((dat$p)^3)/27)

for (i in 1:hertz) {
  if (dat$D[i]>0)
    {
        one1<-(dat$D[i])^0.5
        one2<–0.5*dat$q[i]
        dat$x1[i]<-Math.cbrt(one2+one1)
        dat$x2[i]<-Math.cbrt(one2-one1)
        dat$a[i]<-dat$x1[i]+dat$x2[i]
        rm(one1)
        rm(one2)
    }
 else {
   if (dat$D[i]<0)
   {
     one1<-(dat$D[i])^0.5
     one2<–0.5*dat$q[i]
     dat$x1[i]<-Math.cbrt(one1+one2)
     dat$x2[i]<-Math.cbrt(one1-one2)
     rm(one1)
     rm(one2)
     
     part1<-((-1)*dat$p[i]/3)^3 
    b<-NISTradianTOdeg(acos((-0.5*dat$q[i])/((part1)^0.5)))
    
    a1=2*((-dat$p[i]/3)^0.5)*(cos(NISTdegTOradian(b/3)))
    a2=2*((-dat$p[i]/3)^0.5)*(cos(NISTdegTOradian(120+b/3))) 
    a3=2*((-dat$p[i]/3)^0.5)*(cos(NISTdegTOradian(240+b/3)))
    c<-max(a1, a2, a3)
    dat$a[i]<-c
    rm(c, a1, a2, a3, part1)
    
   }
 } 
}

myvars <- names(dat) %in% c(“x1”, “x2”) 
dat <- dat[!myvars]

##### using polynomial solver for a ####
require(polynom)

for (i in 1:hertz){
  
  p<-polynomial(coef = c(dat$q[i],dat$p[i], 0, 1))
  pz <- solve(p)
  length(pz)
  for (j in 1:length(pz))
       {
    assign(paste(“sol”, j, sep = “”), pz[j]) }
      dat$sol1[i]<-sol1
  dat$sol2[i]<-sol2
  dat$sol3[i]<-sol3
    rm(pz)
    rm(sol1, sol2, sol3)
}


for(i in 1:hertz)
{
z1<-dat$sol1[i]
z2<-dat$sol2[i]
z3<-dat$sol3[i]

dat$realsol1[i]<-Re(z1)
dat$Imsol1[i]<-Im(z1)
dat$realsol1[i]
dat$Imsol1[i]


dat$realsol2[i]<-Re(z2)
dat$Imsol2[i]<-Im(z2)
dat$realsol2[i]
dat$Imsol2[i]

dat$realsol3[i]<-Re(z3)
dat$Imsol3[i]<-Im(z3)
dat$realsol3[i]
dat$Imsol3[i]

max(dat$realsol1[i],dat$realsol2[i],dat$realsol3[i])
dat$valA[i]<-max(dat$realsol1[i],dat$realsol2[i],dat$realsol3[i])
dat$max3lagtime[i]<-dat$moment3[i]/dat$col2[i]
}

myvars <- names(dat) %in% c(“col1”, “col2”,
                            “moment2”, “moment3”,
                            “moment5”,“p”, “q”, “a”,“valA”,  “max3lagtime” ) 
summaryTable <- dat[myvars]

summaryTable ##    col1 col2   moment2    moment3    moment5          p           q
## 1     1 0.05  39.97585   227.4659   262452.8  -754.0533   -2274.659
## 2     2 0.10 115.63242  1134.9434  3160154.1 -1628.0914  -11349.434
## 3     3 0.15 189.03933  2265.2436  8969219.4 -2069.1011  -22652.436
## 4     4 0.20 254.38826  3320.6327 15645921.7 -2167.8465  -33206.327
## 5     5 0.25 312.56747  4242.6024 22015534.6 -2063.4833  -42426.024
## 6     6 0.30 364.54446  5068.9708 27896667.6 -1857.9740  -50689.708
## 7     7 0.35 410.76632  5884.4135 34475044.8 -1751.0421  -58844.135
## 8     8 0.40 452.93659  6660.0190 41266600.0 -1666.8026  -66600.190
## 9     9 0.45 492.31085  7343.4792 47495349.9 -1544.5819  -73434.792
## 10   10 0.50 529.03029  7939.6602 53081193.5 -1395.2721  -79396.602
## 11   11 0.55 562.93084  8509.7643 59150608.6 -1321.6019  -85097.643
## 12   12 0.60 593.42319  9046.5529 66051182.5 -1367.0223  -90465.529
## 13   13 0.65 620.52924  9426.8395 71306547.4 -1358.9127  -94268.395
## 14   14 0.70 644.90293  9647.2450 73292986.2 -1148.2677  -96472.450
## 15   15 0.75 667.44272  9846.4826 73991999.9  -840.1344  -98464.826
## 16   16 0.80 689.19924 10113.5842 75160476.5  -539.6436 -101135.842
## 17   17 0.85 710.79955 10484.7886 78308090.9  -360.7380 -104847.886
## 18   18 0.90 732.54936 10959.1331 83838648.3  -324.6232 -109591.331
## 19   19 0.95 754.19575 11389.8848 89552992.0  -320.5445 -113898.848
## 20   20 1.00 775.68790 11665.0130 94068144.4  -307.2478 -116650.130
##           a     valA max3lagtime
## 1  28.85951 28.85951    4549.317
## 2  43.46502 43.46502   11349.434
## 3  50.20279 50.20279   15101.624
## 4  52.87582 52.87582   16603.163
## 5  53.45273 53.45273   16970.410
## 6  53.04339 53.04339   16896.569
## 7  53.41124 53.41124   16812.610
## 8  53.87871 53.87871   16650.048
## 9  53.91351 53.91351   16318.843
## 10 53.62664 53.62664   15879.320
## 11 53.86497 53.86497   15472.299
## 12 54.90598 54.90598   15077.588
## 13 55.33885 55.33885   14502.830
## 14 54.13317 54.13317   13781.779
## 15 52.21134 52.21134   13128.643
## 16 50.44371 50.44371   12641.980
## 17 49.70186 49.70186   12335.045
## 18 50.11435 50.11435   12176.815
## 19 50.67653 50.67653   11989.352
## 20 50.95577 50.95577   11665.013 ### picking the ramp characteristics based on the rato of third moment and lag time.

shortened_sol<-dat[which(dat$max3lagtime==max(dat$max3lagtime)),]
shortened_sol ##   col1 col2     col3      col4      col5  moment2  moment3  moment5
## 5    5 0.25 312.5675 -4242.602 -22015535 312.5675 4242.602 22015535
##           p         q         D        a               sol1
## 5 -2063.483 -42426.02 124575728 53.45273 -26.72637-8.91137i
##                 sol2        sol3  realsol1    Imsol1  realsol2   Imsol2
## 5 -26.72637+8.91137i 53.45273+0i -26.72637 -8.911368 -26.72637 8.911368
##   realsol3 Imsol3     valA max3lagtime
## 5 53.45273      0 53.45273    16970.41 amplitude<-shortened_sol$a
tau<-((-1*shortened_sol$col2)*(shortened_sol$valA)^3)/shortened_sol$col4
r<-shortened_sol$col2

solution<-data.frame(“Parameters” = c(“r”, “amplitude”, “tau”), 
                     “Values”= c(r, amplitude, tau), 
                     “Units”= c(“seconds”,“mmol m-3”,“seconds” ))
solution ##   Parameters    Values    Units
## 1          r  0.250000  seconds
## 2  amplitude 53.452730 mmol m-3
## 3        tau  8.999479  seconds

Acknowledgements

The author would like to acknowledge Andy Suyker for providing the test dataset. The authors also acknowledges the contributions of Kosana Suvočarev and Tala Awada.


References

Castellvi, F., Martinez-Cob, A., & Perez-Coveta, O. (2006). Estimating sensible and latent heat fluxes over rice using surface renewal. Agricultural and forest meteorology139(1), 164-169.

Castellvi, F., & Snyder, R. L. (2009). Combining the dissipation method and surface renewal analysis to estimate scalar fluxes from the time traces over rangeland grass near Ione (California). Hydrological processes23(6), 842-857.

Castellví, F., & Snyder, R. L. (2009). On the performance of surface renewal analysis to estimate sensible heat flux over two growing rice fields under the influence of regional advection. Journal of hydrology375(3), 546-553.

Castellvi, F., Snyder, R. L., & Baldocchi, D. D. (2008). Surface energy-balance closure over rangeland grass using the eddy covariance method and surface renewal analysis. Agricultural and forest meteorology148(6), 1147-1160.

Katul, G., Hsieh, C. I., Oren, R., Ellsworth, D., & Phillips, N. (1996). Latent and sensible heat flux predictions from a uniform pine forest using surface renewal and flux variance methods. Boundary-Layer Meteorology80(3), 249-282.

Mengistu, M. G., & Savage, M. J. (2010). Surface renewal method for estimating sensible heat flux. Water SA, 36(1), 9-18.

Paw U, K. T., Snyder, R. L., Spano, D., & Su, H. B. (2005). Surface renewal estimates of scalar exchange. Agronomy, 47, 455.

Snyder, R. L., Spano D., Duce K.T., Paw U, Anderson F. E. and Falk M. (2007). Surface Renewal Manual. University of California, Davis, California.

Spano, D., Snyder, R. L., & Duce, P. (2000). Estimating sensible and latent heat flux densities from grapevine canopies using surface renewal. Agricultural and Forest Meteorology104(3), 171-183.

Suvočarev, K., Shapland, T. M., Snyder, R. L., & Martínez-Cob, A. (2014). Surface renewal performance to independently estimate sensible and latent heat fluxes in heterogeneous crop surfaces. Journal of hydrology509, 83-93.

To leave a comment for the author, please follow the link and comment on their blog: R-posts.com.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)