Bulk Aerodynamic Coefficients

Roughness Height for Momentum

Sensible Heat Flux

Wind Speed Range

Aerodynamic Resistances


The transfer of heat, moisture, and momentum (wind) between the atmosphere and the Earth’s surface depends is driven by differences in temperature (for sensible heat flux), water vapor (for latent heat flux), and wind speed (for momentum flux). By analogy to Ohm’s Law for electrical current, the ratio of these differences to the respective fluxes is defined as the aerodynamic resistance:

This toy model implements the method described in sections 14.5 – 14.7 in the Bonan textbook (pages 208-212) to compute the aerodynamic resistance to heat transport as a function of different average wind speeds. The reference height is 10 meters above the surface and the displacement height (\(d\)) is assumed to be zero. The user specifies the roughness height (\(z_{0M}\)) for momentum and the sensible heat flux (\(H\)) on the sliders to the left. The roughness lengths for heat flux (\(z_{0H}\), \(z_{0W}\)) are assumed to be 10% of the specified \(z_{0M}\). For details please see the textbook and/or the lecture notes.

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, using a web programming package called shiny. 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 four important components:

  • A script (resistance.R) that calculates aerodynamic resistances
  • A script (plot.resistance.R) that plots resistance vs wind speed
  • 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!

Actual calculation of aerodynamic resistance


resistance <- function(roughness=1, disp=0, theta=298, heat.flux=0, wind=5){

  # Calculate resistance (transfer coefficients) for turbulent fluxes 
  # using Monin-Obukhov Similarity and Bulk Transfer equations from Bonan 14.6 and 14.7

  # Reference height 
  ref.height <- 10 # meters

  # Some physical constants
  vonKarman <- 0.4
  gravity <- 9.81 # m/s2
  pressure <- 1e5  # Pascals
  C.p <- 1004 # J/K/kg
  R.gas <- 287. # J/K/kg

  # Ideal Gas Law
  rho <- pressure / (R.gas * theta)  # density in kg/m3

  # Nondimensional heights from roughness lengths
  # Roughness lengths for heat & water ~ 10% of z.0
  # (See note in Bonan p 212 following eqn 14.34)
  log.zm <- log((ref.height-disp)/roughness)
  log.zh <- log((ref.height-disp)/(roughness/10))
  log.zw <- log.zh

  # Friction velocity (from Bonan 14.16)
  u.star <- vonKarman * wind / log.zm

  # Obukhov Length (Bonan 14.20)
  Obukhov.L <- -u.star^3 / (vonKarman*gravity/theta * heat.flux / (rho*C.p))

  # Stability parameter (Bonan 14.19)
  zeta <- (ref.height-disp) / Obukhov.L 

  # Similarity functions (Bonan 14.26 and 14.27)
  if (zeta < 0) {              # unstable case
    x <- (1 - 16*zeta)^0.25
    psi.m <- 2 * log((1+x)/2) + log((1+x^2)/2) - 2*atan(x) + pi/2
    psi.h <- 2 * log((1+x^2)/2)
    psi.w <- psi.h
  } else {                    # stable case
    psi.m <- -5 * zeta
    psi.h <- psi.m
    psi.w <- psi.m

  # Resistances from Bonan 14.33 
  r.m <- (log.zm-psi.m)^2 / (vonKarman^2 * wind)
  r.h <- (log.zm-psi.m) * (log.zh-psi.h) / (vonKarman^2 * wind)
  r.w <- (log.zm-psi.m) * (log.zw-psi.w) / (vonKarman^2 * wind)

  return(list(r.m=r.m, r.h=r.h, r.w=r.w))


Program to plot resistance vs mean wind speed


plot.resistance <- function(wind.min=0.1, wind.max=10,
                            roughness=1, disp=0, theta=298, heat.flux=0){

  # Plot aerodynamic resistances for heat and momentum
  # Based on Bonan sections 14.5-14.7 (pp 208-212)

  # Set up arrays for plotting
  wind <- rep(NA,n.steps)
  r.m <- rep(NA,n.steps)
  r.h <- rep(NA,n.steps)

  # Loop through 100 different wind speeds, calculate resistances for each one
  for (i in 1:n.steps){
    wind[i] <- wind.min + (wind.max - wind.min) * i / n.steps
    out <- resistance(wind=wind[i],  
                      heat.flux=heat.flux, roughness=roughness, disp=disp, theta=theta)
    r.m[i] <- out$r.m    
    r.h[i] <- out$r.h    

  # Find max and min resistance for plotting
  plot.max <- max(c(r.m, r.h))
  plot.min <- min(c(r.m, r.h))

  # Draw the plot of resistance vs wind speed
  plot(wind, r.m, typ='l', col='blue', lwd=3, ylim=c(plot.min, plot.max),
       main='Surface Layer Resistances', xlab='Wind Speed at 10 m', ylab='s/m')
  lines(wind, r.h, col='red', lwd=3)

  # Add a nice legend and some text
  legend('topright', c('Momentum', 'Heat & Water'), col=c('blue','red'), lwd=c(3,3))
  text( 0.5*(wind.min+wind.max), 0.7*plot.max, 
        paste('Sensible Heat Flux =',heat.flux, 'W/m2'))


This is the user interface that actually controls this website



# Define user interface for aerodynamic model

  #  Application title
  headerPanel("Bulk Aerodynamic Coefficients"),

  # Sidebar on the left with controls for the user to specify atmospheric properties


    # Set model parameters using sliders
    h4('Roughness Height for Momentum'),            
    sliderInput("roughness", "(z0.m, meters)", 
                min = 0.01, max = 5, value = 1),
    h4('Sensible Heat Flux'),
    sliderInput("heat.flux", "Watts per sq meter", 
                min = -100, max = 100, value = 0, step= 5),
    h4('Wind Speed Range'),
    sliderInput("wind.limits", "meters per second",
                min = 0.1, max = 20, value = c(0.5,5), step= 0.1)

  # Show a tabset in the main panel of the browser that displays model output

      # First tab shows plots of assimilation rate vs light and CO2
               h4('Aerodynamic Resistances'), 

      # Second tab displays a brief model description
      tabPanel("Model Description", 

      # Third tab displays website code
      tabPanel("Website code", 

The “server” that calls other scripts when user does something



# Source required R scripts

shinyServer(function(input, output) {
  #  plot of resistances vs wind
  output$resistance.plot <- renderPlot(
    plot.resistance(roughness=input$roughness, heat.flux=input$heat.flux,