Simulating Surface Energy Budget and Evapotranspiration with RMEP
Yong Yang1, Huaiwei Sun2, Kaixun Wang3
version 3.0, 26-Nov-2021
Introduction
This vignette aims to introduce how to conveniently implement the maximum entropy production(MEP) model with RMEP package. The MEP model (developed by Wang et al., 2009) can calculate the surface heat fluxes with fewer inputs and ensure the surface energy balance compared to classical ET models. The RMEP package is designed to calculate actual evapotranspiration (ET), potential ET, sensible heat flux and ground heat flux, and realize the data pre-processing, statistical analysis of result as well as visualization operations at the same time. The RMEP package supports multiple data formats input in text or netCDF format. The calculating functions including: RMEP, RMEPPET, RMEP_nc, Shum and SShum. The plot functions including: RMEP_Plot_line, RMEP_plot_sca and RMEP_plot_spa. The RMEP_mice function is used for data pre-processing, and the RMEP_sta is used for statistical analysis.
Setting up the environment
Loading the RMEP package first,which contains data and functions used in this analysis. If the RMEP has not been installed, install RMEP package from the github repository.
library(RMEP)Then we show the usage and description of each function below, and provide examples for specific demonstration.
Calculate evapotranspiration(ET) and heat fluxes
Three variables are required at least in theory for the maximum entropy production (MEP) model inputs, including the net radiation (Rn), specific humidity (qs) and surface temperature (Ts). Net long wave radiation (RnL) variable is also introducing in the RMEP for the case of Water-snow-ice surface in the MEP formulations. When calculate potential evapotranspiration (PET), only two variables (Rn and Ts) are needed.
Usage
RMEP(Rn,RnL,qs,Ts,type,I,z)
- This function give a list of ET and heat fluxes
Example
The RMEP compute with inputs of simple vector:
RMEP(Rn=300,RnL=100,qs=0.004,Ts=25,I=800,z=2.5,type=1)
#> $ETMEP
#> [,1]
#> [1,] 2.311116
#>
#> $HMEP
#> [,1]
#> [1,] 123.2513
#>
#> $EMEP
#> [,1]
#> [1,] 65.53513
#>
#> $GMEP
#> [,1]
#> [1,] 111.2136
#>
#> $ef_id
#> [1] 1- After calculating, the actual ET is 2.31 (mm/day), the sensible heat (H), latent heat (E) and grounf heat (G) are 123.25, 65.53, 111.21 (W/m2) respectively, and the effective data number is 1. Parameter “Type” represents the condition of surface land cover, type 1 for bare soil surface or short canopy, 2 for dense canopy and 3 for Water-snow-ice surface, and the default is 1.Parameter “I” (constant) is the media thermal inertial (units:TIU),default is 800 TIU. Parameter “z” (constant) is the theoretical height above surface (unit:m), and default is 2.5 m.
Usage
RMEPPET (Rn,RnL,qs,Ts,type,I,z)
- This function give a list of PET and heat fluxes.
RMEPPET(Rn=300,RnL=100,Ts=25,type=2)
#> [1] "RMEPPET completed!"
#> $PETMEP
#> [,1]
#> [1,] 7.400631
#>
#> $HMEP
#> [,1]
#> [1,] 90.14411
#>
#> $EMEP
#> [,1]
#> [1,] 209.8559
#>
#> $GMEP
#> [,1]
#> [1,] 0
#>
#> $ef_id
#> [1] 1- The PET is 7.386154 mm/day.
Calculate the specific humidity
Usage
Shum (TA,PA,RH)
- This function calculate the actual specific humidity with three variables including: Air temperature(TA, unit: deg C), Relative humidity(RH, unit: percent), and Atmospheric pressure(PA, unit: kPa).
Example
Shum(TA=20,PA=101,RH=50)
#> [1] 0.007311296- The qs is 0.007311296 kg/kg.
Calculate the saturated specific humidity
Usage
SShum (TS)
- This function calculate the saturated specific humidity using surface temperature(Ts, unit: deg C).
Example
SShum(Ts=20)
#> [1] 0.01464032- The saturated specific humidity is 0.01464032 kg/kg.
Calculate the actual ET with Netcdf format input
Usage
RMEP_nc(Rn,RnL,qs,Ts,I,z,type)
- This function calculate actual ET (Unit:W/m2) with Netcdf format input. Empirical parameter I is a function of soil moisture, parameter z is the distance above the surface and can be obtained according to land cover types. We use the monthly Global Land Data Assimilation System (GLDAS) data of January, 2000 for demonstration.
Example
## Before calculating, load the variables with Netcdf format input
## We use the internal data of RMEP package for example
Rn <- Rn;RnL <- RnL;qs <- qs;Ts <- Ts; ## load the necessary inputs
sm <- sm ## soil moisture data
z <- z ## parameter z is provided in the internal RMEP package
I <- sqrt(800^2+sm/100*1557^2) ## calculate the thermal inertia (I)
type <- ifelse(z>4,2,1) ## define type
ETdemo <- RMEP_nc(Rn, RnL, qs, Ts, I, z, type) ## ET calculation
## Show the spatial distribution of the calculated ET using the "RMEP_plot_spa" function
## The RMEP_plot_spa will be introduced later.
lon <- lon;lat <- lat; ## longitude and lattitude
RMEP_plot_spa(lon,lat,ETdemo,main = "MEP ET in 2000/01") 
- The spatial distribution of global land ET in January, 2000 is shown as the map above.
Data pre-processing with the input data
Usage
RMEP_mice(data)
- This function deal with the missing data of the inputs for the maximum entropy production(MEP) model.Input the data with missing value and automatically interpolate it into complete data. The function creates multiple imputations (replacement values) for multivariate missing data. There are six univariate imputation methods available including:Predictive mean matching(“pmm”), Weighted predictive mean matching(“midastouch”), Classification and regression trees (“cart”), Unconditional mean imputation (“mean”), Imputation of quadratic terms(“quadratic”) and Linear regression, predicted values (“norm.predict”). Users can select method for imputations and the default is Predictive mean matching(“pmm”).
Example
## We use the built-in dataset "airquality" in R for example
airquality ## view the data with missing value
RMEP_mice(airquality) ## imputation with the default predictive mean matching ("pmm") method- The processed data “airquality” is complete without missing value.
Statistical analysis
Usage
RMEP_sta(MEP,obs)
- This function implements statistical analysis between MEP predicted and observed values, providing the goodness-of-fit statistics to evaluate the results. Quantitative statistics included are: Mean Error (ME), Mean Absolute Error (MAE), Root Mean Square Error (RMSE), Normalized Root Mean Square Error (NRMSE), Pearson product-moment correlation coefficient (r), Spearman Correlation coefficient (r.Spearman), Coefficient of Determination (R2), Ratio of Standard Deviations (rSD), Nash-Sutcliffe efficiency (NSE), Modified Nash-Sutcliffe efficiency (mNSE), Relative Nash-Sutcliffe efficiency (rNSE), Index of Agreement (d), Modified Index of Agreement (md), Relative Index of Agreement (rd), Coefficient of Persistence (cp), Percent Bias (pbias), Kling-Gupta efficiency (KGE), the coef. of determination multiplied by the slope of the linear regression between ‘sim’ and ‘obs’ (bR2), and volumetric efficiency (VE).
Example
## We use the built-in dataset "datademo" in RMEP package for example.
## The "datademo" data contains 8 variables during 1-30, June 2016 with half-hourly time step.
## Load example data first
data("datademo")
head(datademo) ## view the first 6 rows of data ## Calculate the latent heat with the example data
Data <- datademo
qs <- Shum(TA <- Data$TA,PA <- Data$PA,RH <- Data$RH) ## Calculate specific humidity
Rn <- Data$NETRAD;RnL <- Data$NETRAD ## RnL is not used in type=1 so let it equals to Rn
Ts <- Data$TS ## Surface temperature
output <- RMEP(Rn,RnL,qs,Ts,type=1) #The calculated results are obtained in the "output" list
## All the calculated variables can be analyzed, including ET and heat fluxes.
## The MEP calculated latent heat is evaluted foe example
MEP <- t(output$EMEP)[,1] ## MEP calculated latent heat flux
obs <- Data$LE ## Observations
RMEP_sta(MEP,obs)
#> [,1]
#> ME -24.28
#> MAE 29.64
#> MSE 1618.49
#> RMSE 40.23
#> NRMSE % 37.80
#> PBIAS % -25.70
#> RSR 0.38
#> rSD 0.92
#> NSE 0.86
#> mNSE 0.67
#> rNSE -132.56
#> d 0.96
#> md 0.83
#> rd -35.21
#> cp -0.46
#> r 0.95
#> R2 0.91
#> bR2 0.75
#> KGE 0.73
#> VE 0.69- The model performance of MEP estimate is given. For example, the Root Mean Square Error (RMSE) is 40.23, the Nash-Sutcliffe efficiency (NSE) is 0.86 and the Coefficient of Determination (R2) is 0.91.
Visualization and plot
The RMEP package provides plot functions including: line plot (RMEP_plot_line), scatter plot (RMEP_plot_sca) and spatial map (RMEP_plot_spa).
Usage
RMEP_plot_line(timeseries)
- This function draw a line plot between the MEP estimations and observations. The input data for this function contains three vectors named: time, MEP, obs, which represents time series, MEP estimations and observations, respectively.
Example
## We use the MEP and obs data obtained from the datademo in the previous example
time <- 1:length(Data$year)
timeseries <- data.frame(time,MEP,obs)
RMEP_plot_line(timeseries) 
- The blue line represents the MEP estimations and red line represents observations.
Usage
RMEP_plot_sca(obs,MEP)
- This function draw a scatter plot between the MEP estimations and observations. The coefficient of determination (R2), the regression equation, the 1:1 line (black dash) and the regression line (red solid) are shown in the plot.
Example
RMEP_plot_sca(obs,MEP)
- The R2 is 0.91 and the regression equation is y=0.882x-13.1.
Usage
RMEP_plot_spa(lon,lat,MEP,main,legends)
- This function draw a spatial map of results. It can draw a spatial map of ET, energy fluxes as well as any variables. The longitude and latitude information are required.
Example
An example was given in the front using the “RMEP_nc” function with ETdemo data. At the same time, we put a built-in data of global land mean annual ET during years 1978-2018 in this package. The annual ET is shown below for example.
data("ETannual") ## load the data
annual <- ETannual ## global mean annual ET during 1978-2018
RMEP_plot_spa(lon,lat,annual,main="Global mean annual ET",legends="mm/yr")
The parameter z mentioned above is also provided in RMEP package for reference,as shown in the plot below.
RMEP_plot_spa(lon,lat,z,main="Z value for MEP",legends="meter")
Summary
As shown in the above, the RMEP package can basically realize the whole process from calculation to analysis, including data pre-processing, statistical evaluation and visualization. We add some built-in data for demonstration, including the “datademo”, “z” value as well as other inputs for MEP model. The global land mean annual ET data is also available in the built-in data within the RMEP package. We hope that the developed RMEP package can help users to implement the MEP model conveniently.
Software details
This tutorial was built under:
#> [1] "x86_64-w64-mingw32/x64 (64-bit)"
#> [1] "R version 4.1.2 (2021-11-01)"
#> [1] "RMEP 2.0"References
[1] Wang,J., & Bras,R.L..(2009). A model of surface heat fluxes based on the theory of maximum entropy production. Water Resources Research, 45(11), 130-142.
[2] Yang, Y., Sun, H., and Wang, J. An R tool for Capturing Dynamics of Actual Evapotranspiration with MEP model and its application in Amazon, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-1034, https://doi.org/10.5194/egusphere-egu2020-1034, 2019.
Version history
- v0.3: Nov 2021
- v0.2: Mar 2021
- v0.1: Dec 2020