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.
This website is controlled using the R package “shiny.” There are four important components:
Scroll down or click links in the list above to read all about it!
resistance.R:
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))
}
plot.resistance.R:
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
n.steps=100
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'))
}
ui.R:
library(shiny)
library(markdown)
library(knitr)
# Define user interface for aerodynamic model
shinyUI(pageWithSidebar(
# Application title
headerPanel("Bulk Aerodynamic Coefficients"),
# Sidebar on the left with controls for the user to specify atmospheric properties
sidebarPanel(
img(src='eddies.png'),
# 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
mainPanel(
tabsetPanel(
# First tab shows plots of assimilation rate vs light and CO2
tabPanel("Output",
h4('Aerodynamic Resistances'),
plotOutput('resistance.plot')),
# Second tab displays a brief model description
tabPanel("Model Description",
includeHTML('doc/aerodynamic.description.html')),
# Third tab displays website code
tabPanel("Website code",
includeMarkdown('doc/website.code.md'))
)
)
))
server.R:
library(shiny)
# Source required R scripts
source('model/resistance.R')
source('model/plot.resistance.R')
shinyServer(function(input, output) {
# plot of resistances vs wind
output$resistance.plot <- renderPlot(
plot.resistance(roughness=input$roughness, heat.flux=input$heat.flux,
wind.min=input$wind.limits[1],
wind.max=input$wind.limits[2]))
})