**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 H_{2}O 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 (CH_{4}) and carbon dioxide (CO_{2}) 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.

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, CO_{2} and CH_{4}.

**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_{i }-V_{i-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 2^{nd}, 3^{rd} and 5^{th} 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 H_{2}O concentrations from g per m^{3} to mmol per m^{3}.

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 *q*. *A* (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 meteorology*, *139*(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 processes*, *23*(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 hydrology*, *375*(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. A*gricultural and forest meteorology*, *148*(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 Meteorology*, *80*(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 Meteorology*, *104*(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 hydrology*, *509*, 83-93.

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