Soil Hydraulic Properties

Hydraulic Properties vs. Water Content

Properties of Specified Soil

MODEL DESCRIPTION

This is a very simple model of soil hydraulic properties intended to be used for teaching and learning about land-surface climate.

Soil matric potential (“suction”) and hydraulic conductivity depend strongly on mopisture content as described by Clapp and Hornberger (1978) and on soil texture as described by Cosby et al (1984). This website allows the user to explore those relationsips.

Soil texture is characterized as a percentage content of sand, silt, and clay. Cosby et al (1984) related porosity (\(\theta_{sat}\)), hydraulic conductivity (\(K_{sat}\)) and matric potential at saturation (\(\Psi_{sat}\)) to soil texture. This calculator usses the values rom their Table 4. Following Bonan’s textbook, we find the volumetric water content at field capacity and wilt point by setting soil suction to 1000 mm for field capacity and 150,000 mm for wilt point.

Capp and Hornberger (1978) used power laws to estimate hydraulic conductivity and suction from the saturated values given volumetric water content. Given \(\theta < \theta_{sat}\), \[ \begin{aligned} \Psi = \Psi_{sat} (\frac{\theta}{\theta_{sat}})^{-b} \end{aligned} \] \[ \begin{aligned} k = k_{sat} (\frac{\theta}{\theta_{sat}})^{2b+3} \end{aligned} \]

The program that does the calculation is very simple, and the code that controls this website is surprisingly simple too! It is all written in the programming language R. You can read all about it on the “Website Code” tab to the right.

How this website works (including all the code!)

This website is controlled using the R package “shiny.” There are six important components:

  • A module that calculates soil parameters by texture (cosby.R)
  • A module that applies soil properties at a given moisture (clapp.and.hornberger.R)
  • A module that plots hydraulic properties vs soil moisture (plot.soil.R)
  • A module that crates a formatted table of soil properties (report.output.R)
  • The user interface (ui.R) that controls the sliders and passes inputs to the model
  • The “server” (server.R) that responds to user interface events in the web browser

Scroll down or click links in the list above to read all about it!

This is the main program that calculates soil properties for a given texture

cosby.R:

cosby <- function(percent.clay=18, percent.sand=43) {

  # Calculate soil hydraulic parameters based on Cosby et al (1984)
  # (Coefficients are derived from Cosby et al Table 4)

  # vwc is volumetric water content (m^3 of water per m^3 of soil)
  # suction is matric.potential (soil suction)

  # Density of water
  rho.water <- 1000. # kg per m^3

  # Acceleration of gravity
  gravity <- 9.81 # meters per second

  # Make sure texture percentages sum to 100
  percent.silt <- 100. - (percent.clay + percent.sand) 

  ##################

  # Coefficients from Cosby et al (1984) Table 4:
  bee.intercept <- 3.10
  d.bee.d.clay <- 0.1570
  d.bee.d.sand <- -0.003

  suction.sat.intercept <- 1.54
  dLog.suction.sat.d.sand <- -0.0095
  dLog.suction.sat.d.silt <- 0.0063

  k.sat.intercept <- -0.60
  dLog.k.sat.d.sand <- 0.0126
  dLog.k.sat.d.clay <- -0.0064

  porosity.intercept <- 50.5
  d.porosity.d.clay <- -0.037
  d.porosity.d.sand <- -0.142

  ##################

  # Clapp and Hornberger "b" exponent
  bee = bee.intercept + d.bee.d.clay * percent.clay + d.bee.d.sand * percent.sand

  # Saturated soil water suction psi_s (mm)
  suction.sat = 10 * (10^(suction.sat.intercept + 
                          dLog.suction.sat.d.sand * percent.sand +
                          dLog.suction.sat.d.silt * percent.silt))

  # Saturated soil hydraulic conductivity K_s (mm per hour)
  hydraulic.cond.sat = 25.4 * 10^(k.sat.intercept +
                                    dLog.k.sat.d.sand * percent.sand +
                                    dLog.k.sat.d.clay * percent.clay)

  # Volumetric soil water concentration at saturation point θs (porosity)
  vwc.sat = 0.01 * (porosity.intercept + 
                      d.porosity.d.clay * percent.clay +
                      d.porosity.d.sand * percent.sand)

  # Volumetric soil moisture concentration at Field Capacity (Bonan Fig 9.8)
  suction.field.capacity <- 1000.   # mm
  vwc.field.capacity <- vwc.sat * (suction.field.capacity / suction.sat) ^ (-1/bee)

  # Volumetric soil moisture concentration at wilting point (Bonan Fig 9.8)
  suction.wilt <- 150000   # mm
  vwc.wilt <- vwc.sat * (suction.wilt / suction.sat) ^ (-1/bee)

  return(list(bee=bee, suction.sat=suction.sat, hydraulic.cond.sat=hydraulic.cond.sat,
                vwc.sat=vwc.sat, vwc.field.capacity=vwc.field.capacity, vwc.wilt=vwc.wilt))

}

This is the program that calculates properties for a specified moisture

clapp.and.hornberger.R:

clapp.and.hornberger <- function(vwc=.4, percent.sand=43, percent.clay=18){ 

  # Use the method of Clapp and Hornberger (1978) to estimate the soil
  # matric potential ("suction") and hydraulic conductivity of a soil 
  # given the texture (%sand and %clay) and the volumetric water content (vwc)

  # First compute soil properties from texture using Cosby et al (1984)
  soil <- cosby(percent.sand=percent.sand, percent.clay=percent.clay)

  # soil suction as a function of vwc
  suction <- soil$suction.sat * (vwc / soil$vwc.sat)^-soil$bee

  # hydraulic conductivity as a function of vwc
  hydraulic.cond <- soil$hydraulic.cond.sat * 
                    (vwc / soil$vwc.sat)^(2 * soil$bee + 3)

  return(list(vwc=vwc, percent.sand=percent.sand,
              percent.silt=(100-(percent.sand+percent.clay)),
              percent.clay=percent.clay,
              suction=suction, hydraulic.cond=hydraulic.cond,
              bee=soil$bee, suction.sat=soil$suction.sat, 
              hydraulic.cond.sat=soil$hydraulic.cond.sat,
              vwc.sat=soil$vwc.sat, vwc.field.capacity=soil$vwc.field.capacity, 
              vwc.wilt=soil$vwc.wilt))

}

This is the program that calculates properties for a specified moisture

report.output.R:

report.output <- function(soil){
  attach(soil)
  output.df <- data.frame(Sand=c(format(percent.sand,digits=3), '%'),
                          Silt=c(format(percent.silt,digits=3), '%'),
                          Clay=c(format(percent.clay,digits=3), '%'),
                          VWC=c(format(vwc,digits=3), 'm3/m3'),
                          Porosity=c(format(100*vwc.sat,digits=3), '%'),
                          Matric.Potential=c(format(suction,digits=4), 'mm'),
                          Hydraulic.Conductivity=c(format(hydraulic.cond,digits=4), 'mm/hr'),
                          b.exponent=c(format(bee,digits=4), ' '),
                          Matric.Potential.at.Saturation=c(format(suction.sat,digits=4),'mm'),
                          Hydraulic.Conductivity.at.Saturation=c(format(hydraulic.cond.sat, digits=4), 'mm/hr'),
                          Field.Capacity=c(format(vwc.field.capacity, digits=3), 'm3/m'),
                          Wilt.Point=c(format(vwc.wilt, digits=3), 'm3/m'))
  detach(soil)
  return(t(output.df))
}

This is the program that plots the soil properties vs moisture

plot.soil.R:

plot.soil <- function(vwc=0.4, percent.clay=27, percent.sand=43){

  # Plot matric potential and hydraulic conductivity as a function of 
  # volumetric percent of saturation

  # Calculate the osil hydraulic parameters based on texture from Cosby et al (1984)
  soil <- cosby(percent.sand=percent.sand, percent.clay=percent.clay)

  # Use Clapp and Hornberger power laws to get hydraulic properties vs water content
  suction <- soil$suction.sat * seq(.05,1,.05) ^ -soil$bee
  cond <- soil$hydraulic.cond.sat * seq(.05,1,.05)^(2 * soil$bee + 3)

  # Make two side-to-side plots of properties vs water content

  # Remember original plot parameters
  orig.par <- par(no.readonly=TRUE)

  # Set up a multifigure page with two rows in one column
  par(mfrow=c(1,2))

  # Plot the matric potential ("suction") as a function of saturation
  plot(seq(5,100,5), suction, log='y', typ='b', xlim=c(0,100),
       cex.axis=1.3, cex.lab=1.2, lwd=3,
       xlab='percent of saturation', ylab='suction (mm)')

  # Place a vertical line at the current percent of saturation
  percent.sat <- 100 * (vwc / soil$vwc.sat)
  abline(v=percent.sat, col='red', lwd=2)

  # Plot the hydraulic conductivity as a function of saturation
  plot(seq(5,100,5), cond, log='y', typ='b', xlim=c(0,100),
       cex.axis=1.3, cex.lab=1.2, lwd=3,
       xlab='percent of saturation', ylab='conductivity (mm/s)')
  abline(v=percent.sat, col='red', lwd=2)

  # Restore original plot parameters
  par(orig.par)

}

This is the user interface that actually controls this website

ui.R:

library(shiny)
library(markdown)

# Define UI for the bucket model of surface hydrology
shinyUI(pageWithSidebar(

  #  Application title
  headerPanel("Soil Hydraulic Properties"),

  # Sidebar with sliders that control various available options
  sidebarPanel(
    img(src='soil.jpg'),

    wellPanel(
      sliderInput("percent.sand", "Percent Sand:", 
                min = 0, max = 100, value = 43, step=1),

      sliderInput("percent.clay", "Percent Clay:", 
                  min = 0, max = 100, value = 27, step=1),

      sliderInput("vwc", "Volumetric Water Content (m3/m3):",
                  min=0, max=0.5, value=0.35, step=0.01)
    )

  ),

   # Main panel consists of a tabset displays model output or model description
   mainPanel(
     tabsetPanel(
       tabPanel("Output", 
                h3("Hydraulic Properties vs. Water Content"),
                plotOutput("hydraulic.properties.plot"),
                h4("Properties of Specified Soil"),
                tableOutput('soil.properties'),
                tags$style(type="text/css", ".tab-content { overflow: visible; }")), 
       tabPanel("Description", 
                includeHTML('doc/soil.description.html')),
       tabPanel("Website Code",
                includeMarkdown('doc/soil.code.md'))
     )
  )
))

Here is the server module that responds to user interface events

server.R:

library(shiny)

# Source required R scripts
source('model/cosby.R')
source('model/clapp.and.hornberger.R')
source('model/plot.soil.R')
source('model/report.output.R')

shinyServer(function(input, output) {

  # Calculate the full suite of hydraulic properties for the soil
  output$soil.properties <- 
    renderTable(expr=report.output(clapp.and.hornberger(
                                   vwc=input$vwc,
                                   percent.sand=input$percent.sand,
                                   percent.clay=input$percent.clay)),
                 include.colnames=FALSE)

  # Plot the values of suction and conductivity vs soil moisture
  output$hydraulic.properties.plot <- 
    renderPlot(
      plot.soil(vwc=input$vwc,
                percent.sand=input$percent.sand,
                percent.clay=input$percent.clay)
    )

})