11 Field Soil Water Balance – R for Spatial Analysis & Visualization

11.1 Field Soil Water Balance

We will develop a 2-layer simple soil water balance modeling for simulating 1-D fluxes in the rootzone using USCRN station at Stillwater, OK, USA as an example. Let us extract sample data from the sample station for 2021.

# Yearly data from the sample station
CRNdat = read.csv(url("https://www.ncei.noaa.gov/pub/data/uscrn/products/daily01/2021/CRND0103-2021-OK_Stillwater_2_W.txt"), header=FALSE,sep="")

# Data headers
headers=read.csv(url("https://www.ncei.noaa.gov/pub/data/uscrn/products/daily01/headers.txt"), header=FALSE,sep="")

# Column names as headers from the text file
colnames(CRNdat)=headers[2,1:ncol(CRNdat)]

# Replace fill values with NA
CRNdat[CRNdat == -9999]=NA
CRNdat[CRNdat == -99]=NA
CRNdat[CRNdat == 999]=NA

# View data sample
library(kableExtra)
dataTable = kbl(head(CRNdat,6),full_width = F)
kable_styling(dataTable,bootstrap_options = c("striped", "hover", "condensed", "responsive"))

WBANNO	LST_DATE	CRX_VN	LONGITUDE	LATITUDE	T_DAILY_MAX	T_DAILY_MIN	T_DAILY_MEAN	T_DAILY_AVG	P_DAILY_CALC	SOLARAD_DAILY	SUR_TEMP_DAILY_TYPE	SUR_TEMP_DAILY_MAX	SUR_TEMP_DAILY_MIN	SUR_TEMP_DAILY_AVG	RH_DAILY_MAX	RH_DAILY_MIN	RH_DAILY_AVG	SOIL_MOISTURE_5_DAILY	SOIL_MOISTURE_10_DAILY	SOIL_MOISTURE_20_DAILY	SOIL_MOISTURE_50_DAILY	SOIL_MOISTURE_100_DAILY	SOIL_TEMP_5_DAILY	SOIL_TEMP_10_DAILY	SOIL_TEMP_20_DAILY	SOIL_TEMP_50_DAILY	SOIL_TEMP_100_DAILY
53926	20210101	2.622	-97.09	36.12	2.1	0.2	1.2	0.8	23.0	3.51	C	1.0	-0.3	0.3	93.1	70.1	86.2	0.474	0.440	0.458	0.450	0.429	3.9	4.2	5.4	8.3	9.9
53926	20210102	2.622	-97.09	36.12	6.9	-2.5	2.2	1.8	0.0	7.02	C	9.9	-4.7	1.1	94.2	57.3	78.7	0.466	0.443	0.459	0.448	0.432	4.6	4.7	5.5	8.0	9.6
53926	20210103	2.622	-97.09	36.12	9.1	-3.7	2.7	1.3	0.0	10.56	C	12.9	-6.2	-0.2	94.2	56.4	83.5	0.436	0.440	0.460	0.449	0.429	4.2	4.7	5.5	7.8	9.4
53926	20210104	2.622	-97.09	36.12	14.7	-3.2	5.7	4.3	0.0	11.13	C	16.8	-6.0	2.0	94.0	46.2	76.4	0.423	0.428	0.461	0.450	0.430	4.6	5.0	5.6	7.6	9.2
53926	20210105	2.622	-97.09	36.12	16.2	-3.6	6.3	6.2	0.0	9.99	C	17.9	-6.5	3.9	94.0	22.0	59.5	0.412	0.417	0.463	0.453	0.431	4.8	5.2	5.7	7.6	9.1
53926	20210106	2.622	-97.09	36.12	9.6	4.8	7.2	7.5	0.3	1.03	C	11.4	4.3	7.1	92.1	45.8	66.7	0.395	0.412	0.459	0.450	0.431	6.2	6.1	6.2	7.6	9.1

11.1.1 Exploratory Data Analysis for the Sample Station

 # Create Date array for the station
 library(zoo)
 library(cowplot)
 library(ggplot2)
  
  # Time (in days)
  datetime=as.Date(as.character(CRNdat$LST_DATE), format= "%Y%m%d")
  
  # Surface soil moisture (in m3/m3)
  SSM=na.approx(CRNdat$SOIL_MOISTURE_5_DAILY,na.rm = FALSE)  
  # Soil moisture at 50 cm depth (in m3/m3)
  SM50=na.approx(CRNdat$SOIL_MOISTURE_50_DAILY,na.rm = FALSE)  

  # Precipitation (in mm)
  PCP=CRNdat$P_DAILY_CALC           
  PCP[PCP<0]=0       # Negative (flagged values are replaced by 0)
  PCP[is.na(PCP)]=0  # NA values are replaced by 0)
  
  # Create a data frame of SM-PCP dataset for plotting
  df=data.frame(Date=datetime, PCP=PCP, SSM=SSM, SM50=SM50)

  # Plot for surface soil moisture 
  p1 = ggplot(df) +
    geom_line(aes(Date, SSM, color = "SSM")) +
    geom_line(aes(Date, SM50, color = "SM50")) +
    scale_y_continuous(position = "left",
                       limits = c(0.05,0.475),
                       expand = c(0,0)) +
    scale_color_manual(values = c("steelblue", "cyan4")) +
    guides(x = guide_axis(angle = 90)) +
     scale_x_date(date_breaks = "month",
                     limits = c(datetime[1], 
                                datetime[365]),
                 expand = c(0,0)) +
    labs(y = "Soil moisture [v/v]",
         x = "Date") +
    theme_minimal() +
    theme(axis.title.y.left = element_text(hjust = 0),
          legend.position = "bottom",
          legend.justification = c(0.1, 0.5),
          legend.title = element_blank())

  # Plot for precipitation 
  p2= ggplot(df) +
    geom_col(aes(Date, PCP, fill = "Total Daily Precipitation")) +
    scale_y_reverse(position = "right",
                    limits = c(100,0),
                    breaks = c(0,50,100),
                    labels = c(0,50,100),
                    expand = c(0,0)) +
    scale_x_date(date_breaks = "month",
                     limits = c(datetime[1], datetime[365]),
                 expand = c(0,0)) +
    scale_fill_manual(values = c("sienna1")) +
    guides(x = guide_axis(angle = 90)) +
    labs(y = "Precipitation [mm]", x = "") +
    theme_minimal() +
    theme(axis.title.y.right = element_text(hjust = 0),
          legend.position = "bottom",
          legend.justification = c(0.75, 0.5),
          legend.title = element_blank())

aligned_plots = cowplot::align_plots(p1, p2, align = "hv", axis = "tblr")
p3=ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])
print(p3)

11.2 Low-pass Filtering of Meteorological Perturbations Through Soil Profile

# Profile depth versus soil temperature
Tempdf=data.frame(T5=na.approx(CRNdat$SOIL_TEMP_5_DAILY,na.rm = FALSE), 
              T10=na.approx(CRNdat$SOIL_TEMP_10_DAILY,na.rm = FALSE), 
              T20=na.approx(CRNdat$SOIL_TEMP_20_DAILY,na.rm = FALSE), 
              T50=na.approx(CRNdat$SOIL_TEMP_50_DAILY,na.rm = FALSE),
              T100=na.approx(CRNdat$SOIL_TEMP_100_DAILY,na.rm = FALSE))

df=rbind(data.frame(Temp=Tempdf$T5, Depth=5,Date= datetime),
data.frame(Temp=Tempdf$T10, Depth=10, Date=datetime),
data.frame(Temp=Tempdf$T20, Depth=20,Date= datetime),
data.frame(Temp=Tempdf$T50, Depth=50,Date=datetime),
data.frame(Temp=Tempdf$T100, Depth=100, Date=datetime))

ggplot(df, aes(x = Date, y = Depth)) +
  geom_contour_filled(aes(z = Temp), bins = 6,   stat = "contour_filled") +
  scale_fill_brewer(palette = "YlOrRd",direction = 1) +
  labs(x = "Date", y = "Depth (cm)", fill = "Temperature (°C)") +
  theme_minimal() + 
  scale_y_reverse()

# Profile depth versus soil moisture
SMpdf=data.frame(SM5=na.approx(CRNdat$SOIL_MOISTURE_5_DAILY,na.rm = FALSE), 
              SM10=na.approx(CRNdat$SOIL_MOISTURE_10_DAILY,na.rm = FALSE), 
              SM20=na.approx(CRNdat$SOIL_MOISTURE_20_DAILY,na.rm = FALSE), 
              SM50=na.approx(CRNdat$SOIL_MOISTURE_50_DAILY,na.rm = FALSE),
              SM100=na.approx(CRNdat$SOIL_MOISTURE_100_DAILY,na.rm = FALSE))

df=rbind(data.frame(SM=SMpdf$SM5, Depth=5,Date= datetime),
data.frame(SM=SMpdf$SM10, Depth=10, Date=datetime),
data.frame(SM=SMpdf$SM20, Depth=20,Date= datetime),
data.frame(SM=SMpdf$SM50, Depth=50,Date=datetime),
data.frame(SM=SMpdf$SM100, Depth=100, Date=datetime))

ggplot(df, aes(x = Date, y = Depth)) +
  geom_contour_filled(aes(z = SM), bins = 6,   stat = "contour_filled") +
  scale_fill_brewer(palette = "YlGnBu") +
  labs(x = "Date", y = "Depth (cm)", fill = "Soil Moisture") +
  theme_minimal() + 
  scale_y_reverse()

Under hydrologic equilibrium, the analytical solution of the 1-dimensional vertical water balance (differential) equation can be approximated using the relationship between the temporal changes in the rootzone soil moisture (RZSM) and the difference in SSM vis-à-vis the antecedent rootzone conditions (Wagner et al., 1999, Albergel et al., 2008, Manfreda et al., 2014). A simple recursive exponential LP filter to simulate the 1-dimensional first-order infinite-impulse response of RZSM to temporal variability in SSM is given as follows:

\[ SM_t^{'}= SM_{t-1}^{'}+ D_p(SSM_t-SM_{t-1}^{'}) \tag{11.1}\]

where, \(t \ge1\); when \(t=1\), \(SM_{t=0}^{'}=SSM_{t=0}\).

\(SM_{t}^{'}\) is the temporally smoothed (filtered) SSM at time \(t\) .

Parameter \(D_p\) is the exponential filter smoothing factor of the soil profile \([-]\), also called the pseudo-diffusivity coefficient; \(D_p\) ∈ [0,1]. It represents the effective influence of various bio-geo-physical controls such as vegetation, topography, hydroclimatology and pedological characteristics (soil profile thickness, effective soil hydraulic characteristics, etc.) on vertical fluxes between the soil surface and the rootzone.

Eq. 1 provides a first-order approximation of the infinite-impulse response of RZSM to temporal variability in SSM as \(SM^{'}\), and yields the exponential weighted average of the antecedent observations using weights proportional to the terms of the geometric progression: \(1, (1-D_p), (1-D_p)^2, ….(1-D_p)^n\), \(-\) the discrete form of an exponential function (Hunter, 1986; Perry, 2010). This method makes two key assumption:

The lateral moisture flux to/from the rootzone is negligible, and,
The groundwater table is significantly deeper than the rootzone depth.

11.2.1 Origins of Low-pass Filter Relationship in Vertical Soil Fluxes

The origins of the low-pass filter effect of soil profile on surface variability in soil moisture can be traced to the original equations for vertical flux of water int he unsaturated soil profile. We recall the Diffusivity version of Richard’s equation as follows:

\[ \frac{\partial\theta} {\partial t}= \frac{\partial}{\partial z}(D(\theta)\frac{\partial\theta}{\partial z})-\frac{\partial K(\theta)}{\partial z} \tag{11.2}\]

where, \(D\) and \(K\) are the hydraulic diffusivity and conductivity of the soil profile, and are a function of the water content of the soil profile, \(z\) is the height above the bottom of the soil profile and \(\theta\) is the soil moisture content. The right hand side of this equation contains two terms, where \(\partial/\partial z (D(\theta)\partial/\partial z)\) accounts for the flux due to the suction gradient in the profile, while \(\partial K(\theta)/\partial z\) pertains to the role of gravity in the vertical flow of water. Given the equivalent suction exerted by the atmosphere is several order of magnitude compare to the gravitational head gradient, we can simplify the preceeding equation as:

\[ \frac{\partial\theta} {\partial t}= \frac{\partial}{\partial z}(D(\theta)\frac{\partial\theta}{\partial z}) \tag{11.3}\]

First integral of this equation can be written as: \[ \frac{d\theta}{dt}=D(\theta)\frac{d\theta}{dz}+C \tag{11.4}\] where, \(C\) is the constant of integration. Since there is no flow at the bottom boundary of the soil profile, \(d\theta/dz=0\) . Hence, \(C=0\). Writing the preceding equation in discrete form, we get:

\[ \frac{\Delta\theta}{\Delta t}=D(\theta)\frac{\Delta\theta}{\Delta z} \tag{11.5}\]

As the soil reaches the second phase of drying, it enters a phase where the rate of soil drydown is slow, and the rate of loss of water from all depths can be assumed to be the same. Simplifying the previous equation for a simple 2-layer representaiton of soil profile, with \(SSM\) and \(SM\) representing the surface and rootzone soil moisture respectively, we can write Eq. 5 as:
\[ SM_t-SM_{t-1}=D(SM) \times (SSM_{t}-SM_{t-1}) \tag{11.6}\]

or,

\[ SM_t=SM_{t-1}+ D(SM) \times (SSM_{t}-SM_{t-1}) \tag{11.7}\]

Note

This equation is analogous to the low-pass filter equation provided earlier. Given the true form of the relationship between the hydraulic diffusivity (or pseudo-diffusivity) with the soil moisture content is not known in the field, often times, the value of \(D\) and \(D_p\) is assumed to be constant for simplicity.

11.2.2 Propagation of Drying Front in the Soil Profile

We will demonstrate the application of low-pass filter to capture the propagation of drying front from the surface to the deeper profiles.

# Function for implementing LP filter on surface soil moisture
lp_filt=function(SSM, Dp){
  # PARAMETERS---
  #   SSM= Surface soil moisture
  #   Dp= Pseudo-diffusivity coefficient
  # RETURNS---
  #   SSMf=filtered surface soil moisture (excluding the spinup period)

  SSMf=c()                           # Empty array to store filtered SSM 
  buffer=0.01*(max(SSM)-min(SSM))    # To decrease sensitivity to noise in sensor readings (1% of SSM range)

  # Apply recursive filter
  SSMf[1]=SSM[1]                    # Initialize the series
  for (t in 2:length(SSM)) 
  SSMf[t]=ifelse(SSM[t] > (SSM[t-1]+buffer),max(SSM[t], SSMf[t-1]),SSMf[t-1]+Dp*(SSM[t]-SSMf[t-1]))
  return(SSMf) 
}

# Filtered surface soil moisture with varying values of Dp
SSM=SMpdf$SM5   # Surface soil moisture
SSMf_1=lp_filt(SSM, Dp= 0.1)
SSMf_05=lp_filt(SSM, Dp= 0.05)
SSMf_02=lp_filt(SSM, Dp= 0.02)

# Plot surface soil moisture with filtered time series
plot(SSM, type="l", col="gray80", lwd=3, ylab="Value", xlab="Time", 
main="Filtered v/s observed SSM for varying Dp")
lines(SSMf_1, col="red")
lines(SSMf_05, col="purple")
lines(SSMf_02, col="blue")

# Normalized filtered surface soil moisture
SSMf_1norm=(SSMf_1-min(SSMf_1))/(max(SSMf_1)-min(SSMf_1))
SSMf_05norm=(SSMf_05-min(SSMf_05))/(max(SSMf_05)-min(SSMf_05))
SSMf_02norm=(SSMf_02-min(SSMf_02))/(max(SSMf_02)-min(SSMf_02))

# Weighted average of soil moisture at multiple depths to estimate RZSM
RZSM=(SMpdf$SM5*(7.5-0)+    # 0- 7.5 cm depth
SMpdf$SM10*(15-7.5)+        # 7.5-15 cm depth
SMpdf$SM20*(35-15)+         # 15-35 cm depth
SMpdf$SM50*(75-35)+         # 35-75 cm depth
SMpdf$SM100*(100-75))/100   # 75-100 cm depth

rzsm_norm=(RZSM-min(RZSM))/(max(RZSM)-min(RZSM))
ssm_norm=(SSM-min(SSM))/(max(SSM)-min(SSM))

plot(rzsm_norm, type="l", col="gray80", lwd=3, ylab="Value", xlab="Time", 
main="Comparison of RZSM with SSM and filtered SSM (normalized to 0-1)")
lines(ssm_norm, col="red")
lines(SSMf_05norm, col="blue")

11.3 Priestley-Taylor Method for Estimating Potential ET (PET)

This formulation is based on physical principles and can be seen as a simplified version of the Penman-Monteith formula for calculating evapotranspiration. The simplification involves removing the aerodynamic terms from the Penman-Monteith equation and using a constant (\(α\)), which has been empirically derived and is estimated to be 1.26 for open bodies of water, but can range from <1 (humid conditions) to almost 2 (arid conditions). This constant tends to be higher in arid regions or areas experiencing significant water stress. The Priestley-Taylor formula can be used for calculating daily evapotranspiration and can also be applied to smaller time intervals (hourly, for example) if the necessary data is available.

\[ \text { PET }=\left(\frac{\alpha}{\lambda_v \rho_{\mathrm{w}}}\right) \times \frac{\Delta}{\Delta+\gamma} \times \text {R} \]

if \(T_{mean}< 0^oC\)

\[ \Delta=0.3405 \times \mathrm{e}^\left(0.0642 \mathrm{~T}_{\text {mean}}\right) \\ \]

if \(T_{mean}> 0^oC\) \[ \Delta=0.3221 \times \mathrm{e}^\left(0.0803 \mathrm{~T}_{\mathrm{mean}}^{0.8876}\right) \]

where,

\(PET\)= Daily potential evapotranspiration (\(m\))

\(λv\) = Latent heat of vaporization (2260 \(kJ/kg\))

\(ρw\) = Density of water (1000 \(kg/m^3\))

\(γ\) = Psychrometric constant ( 4.95x \(10^{-4}\) \(kg/m^3/^oC\))

\(Rad\) = Daily solar radiation incident on a horizontal earth surface (\(kg/m^2\))

\(Δ\) = Slope of the saturation vapor density curve from the psychrometric chart (\(kg/^oC\))

\(T_{mean}\) = Daily mean air temperature (\(^oC\))

\(\alpha\) = Priestley-Taylor correction factor

# Parameters for Priestley-Taylor formula
lambdaV= 2260                  # Latent heat of vaporization 
rhoW= 1000                     # Density of water 
gamma= 4.95*1e-4               # Psychrometric constant
alpha=1.2                      # Priestley-Taylor correction factor

# PET estimation using Priestley-Taylor method 
 PT_PET=function(Tmean,Rad,alpha){
    
    deltaVal=rep(0, length(Tmean))  # Create empty array
    
    # When Tmean<0
    deltaVal[Tmean<0]=0.3405*exp(0.0642*Tmean[Tmean<0])
    # When Tmean>=0
    deltaVal[Tmean>=0]=0.3221*exp(0.0803*((Tmean[Tmean>=0])^0.8876))
    
    # Priestley-Taylor equation 
    PETval=(alpha/(lambdaV*rhoW))*(deltaVal/(deltaVal+gamma))*Rad
    PETval[PETval<0]=0
    return(PETval) #estimate in m
  }

As we can see from the function above, Priestley-Taylor formula requires mean temperature and solar radiation as inputs. So, we will collect the meteorological dataset namely a) Mean daily temperature b) Total solar radiation from the sample station and estimate the value of total daily potential evapotranspiration (in m). We will also store total daily precipitation, which will be used as model forcing.

 # Observed meteorological forcings 
 Rad=CRNdat$SOLARAD_DAILY*1000        # (as 1 MJ= 1000 KJ)
 Tmean=CRNdat$T_DAILY_MEAN
 PCP=CRNdat$P_DAILY_CALC/10           # convert units from mm to cm

 # Linear gap-filling of missing samples and removing flagged samples
 library(zoo)
 Tmean=na.approx(Tmean,na.rm = FALSE)
 Rad=na.approx(Rad,na.rm = FALSE)
 
 PCP[PCP<0]=0       # Negative (flagged values are replaced by 0)
 PCP[is.na(PCP)]=0  # NA values are replaced by 0)

 # Create data frame for met forcings for the station
 datetime=as.Date(as.character(CRNdat$LST_DATE), format= "%Y%m%d")
 forcings=data.frame(Date=datetime,PCP, Rad,Tmean)
  
 # Apply Priestley-Taylor method for PET estimation using the PT_PET function
 PET=PT_PET(Tmean,Rad,alpha)
 
 # Plot estimated PET time series
 ggplot(data=forcings)+
  geom_line(aes(x = datetime, y = PET), color="steelblue") + 
  xlab("Date") +             # Add labels to the axes
  ylab("PET [in m]")+
  theme_bw()