N-Layer Blackbody Atmosphere

Set Model Parameters

Heat transfer by thermal blackbody radiation

MODEL DESCRIPTION

This is a very simple model of vertical energy transfer and temperatures in the Earth system, intended to be used for teaching and learning about weather and climate.

There are two layers of atmosphere and a single layer of surface. The atmosphere is perfectly transparent to solar radiation, so sunshine warms the surface directly without changing the atmosphere. The warm surface then transfers heat to the overlying atmosphere both by radiation of (thermal) infrared radiation and by upward convective heat flux.

The user can use the sliders at left to manipulate the absorptivity/emissivity of each atmospheric layer, the amount of convective heat flux transferred at each level, the intensity of solar radiation, and the albedo of the planet. You can read about how the website works by selecting the tab “Website Code” to the right of this tab in your browser.

The diagram on the right shows the resulting temperature of the surface and each layer of the atmosphere as well as the magnitude of each heat transfer in Watts per square meter. The diagram updates automatically every time the setting are changed.

The model consists of three equations (conservation of energy in each layer) in three unknowns (the three temperatures). You can read more about the underlying physics of this simple energy balance model by viewing the “Physics” tab to the right of this one.

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 four important components:

  • 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
  • A very simple energy balance model (rce.R) that does the calculations
  • A somewhat longer program that draws (make.layer.diagram.R) and updates (update.values.R) the diagram

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


This is the user interface that actually controls this website

ui.R:

library(shiny)
library(markdown)

# Define UI for slider demo application
shinyUI(pageWithSidebar(

  #  Application title
  headerPanel("2-Layer Atmosphere with Solar, Longwave, & Convection"),

  # Sidebar with sliders that demonstrate various available options
  sidebarPanel(
    h4('Set Model Parameters'),            
    sliderInput("eps1", "Emissivity of Lower Atmosphere", 
                min = 0, max = 1, value = 0.5, step= 0.01),

    sliderInput("eps2", "Emissivity of Upper Atmosphere", 
                min = 0, max = 1, value = 0.5, step= 0.01),

    sliderInput("H1", "Convective Heat Flux from Surface to Lower Atmosphere", 
                min = 0, max = 200, value = 0, step=10),

    sliderInput("H2", "Convective Heat Flux from Lower to Upper Atmosphere", 
                min = 0, max = 200, value = 0, step=10),

    sliderInput("solar", "Solar Constant (W/m2)", 
                min = 800, max = 1500, value = 1368, step=10),

    sliderInput("albedo", "Planetary Albedo", 
                min = 0.05, max = 0.9, value = 0.3, step=0.01)
  ),

   # Main panel consists of a tabset displays model output or model description
   mainPanel(
     tabsetPanel(
       tabPanel("Simulation", 
                h4("A simple model of vertical heat transfer in the Earth's climate system"),
                plotOutput("outputValues"),
                tags$style(type="text/css", ".tab-content { overflow: visible; }")), 
       tabPanel("Model Description", 
                includeMarkdown('doc/model.description.md'))
     )
  )
))

Here is the server module

server.R:

library(shiny)

# Source required R scripts
source('model/rce.R')        
source('model/updateValues.R')        

shinyServer(function(input, output) {

  # Table of all the values in modelOutput
  output$outputValues <- renderPlot({
    source('model/make.layer.diagram.R')    
    updateValues(rce(solar=input$solar, albedo=input$albedo, 
                     eps1=input$eps1, eps2=input$eps2, 
                     Hs=input$H1, HL=input$H2))
  }, height=600)

})

This is the radiative-convective model that solves the heat balance of the 2-layer model

rce.R:

rce = function (solar=1367., albedo=0.3, eps1=0.5, eps2=0.5, Hs=50., HL=50.){

  # solar     Solar constant at top of atmosphere (W/m^2)
  # albedo    Shortwave albedo of planet
  # eps1            emissivity of lower atmosphere
  # eps2            emissivity of upper atmosphere
  # Hs              convective heat flux from surface to lower atmosphere (W/m^2)
  # HL              convective heat flux from lower to upper atmosphere (W/m^2)

  sigma <- 5.67e-8  # Stefan-Boltzmann constant (W m^-2 K^-4)

  # Post-albedo solar power to heat the surface
  S <- solar/4*(1-albedo)

  # transfer matrix contains the energy budget coefficients for sfc and each layer
  transfer <- matrix(c(-1. , eps1, (1.-eps1)*eps2, 
                       eps1, -2.*eps1 , eps1*eps2,  
                       (1.-eps1)*eps2 , eps1*eps2 , -2.*eps2),
                       byrow=TRUE, nrow=3)

  # Right-hand-side vector in the energy budget equations      
  rhs <- matrix(c((Hs-S)/sigma , (HL-Hs)/sigma , -HL/sigma))

  # Solve planetary energy budget
  temps4 <- solve(transfer, rhs)  # invert matrix A and multiply by vector from right-hand-side

  # Extract temperatures by taking fourth root of each element of solution vector 
  Ts <- temps4[1]^.25     
  T1 <- temps4[2]^.25
  T2 <- temps4[3]^.25

  # Compute OLR at TOA
  olr <- sigma*(Ts^4*(1-eps1)*(1-eps2) + eps1*T1^4*(1-eps2) + eps2*T2^4)

  # return all the needed information to the calling function (server.R)
  return(list(Tsfc=Ts, Tlower=T1, Tupper=T2, outgoing.longwave=olr, 
              solar=solar, albedo=albedo, solar.absorbed=S,
              emissivity.lower=eps1, emissivity.upper=eps2,
              convec.surface=Hs, convec.lower=HL,
              emission.surface=sigma*Ts^4,
              emission.lower=eps1*sigma*T1^4,
              emission.upper=eps2*sigma*T2^4))
}

This is the code that draws the nice interactive diagram on the main panel

It uses the R package “diagram.”

make.layer.diagram.R:

library(diagram)
# Draw the base diagram

surface.emission.pos <- 0.9
lower.emission.pos <- 0.8
upper.emission.pos <- 0.9
convec.pos <- 0.17
layer.name.pos <- 0.35
solar.pos <- 0.1
box.pos <- 0.5
temp.pos <- 0.6

upper.height=0.84  # midpoint height of the upper box.
lower.height=0.44  # midpoint height of the lower box.
surface.height=0.04  # midpoint height of the surface box.

upper.color = 'dodger blue'
lower.color = 'light blue'
surface.color = 'light green'

# Set plot margins
bottomMargin <- 0.02
leftMargin <- 0.1
topMargin <- 0.1
rightMargin <- 0.1
par(mfrow=c(1,1), mar=c(bottomMargin,leftMargin,topMargin,rightMargin),                     oma=c(0,0,0,0))

radx=0.47   #horizontal radius of the boxes
rady=0.12   #vertical radius of the boxes
lwd=2   #line width of line surrounding the box.
shadow.size=.004    #relative size of shadow.
arrow.thickness <- 5

######################
# Draw the diagram

# Open a new plot window
openplotmat(xlim = c(-0.1, 1.1))

# Draw the incoming radiation from the Sun
straightarrow(c(0.05,1.1),c(solar.pos,surface.height+rady*1.12), 
              lcol='yellow', arr.pos=1,
              endhead=TRUE, lwd=3*arrow.thickness, arr.width=1)
textplain(c(0.05,1.02), lab=expression(over(S[0],4) * (1-alpha) == phantom(0)), cex=1.2)

# Draw three boxes representing upper & lower atmosphere plus surface
shadowbox(box.type='rect', c(box.pos,upper.height), radx, rady=rady, lwd=lwd, 
          shadow.size=shadow.size,  box.col = upper.color)
shadowbox(box.type='rect', c(box.pos,lower.height), radx, rady=rady, lwd=lwd, 
          shadow.size=shadow.size,  box.col = lower.color)
shadowbox(box.type='rect', c(box.pos,surface.height), radx, rady=rady, lwd=lwd, 
          shadow.size=shadow.size,  box.col = surface.color)

# Label the three boxes with their names
textplain(c(layer.name.pos,upper.height),  adj=c(0.5,0.5), 
         lab = c('Upper','Atmosphere'), cex=1.5)
textplain(c(layer.name.pos,lower.height),  adj=c(0.5,0.5), 
         lab = c('Lower', 'Atmosphere'), cex=1.5)
textplain(c(layer.name.pos,surface.height),  adj=c(0.5,0.5), 
          lab = 'Surface', cex=1.5)

# Convective heat flux (thermals) from surface to lower atmosphere
textplain(c(convec.pos,lower.height+rady/3), lab='convec', adj=c(0.5,0.5))
straightarrow(c(convec.pos,lower.height+rady/2), 
              c(convec.pos,upper.height-rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Convective heat flux (thermals) from lower to upper atmosphere
textplain(c(convec.pos,surface.height+rady/3), lab='convec', adj=c(0.5,0.5))
straightarrow(c(convec.pos,surface.height+rady/2), 
              c(convec.pos,lower.height-rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Radiation emitted from upper to lower
textplain(c(upper.emission.pos,upper.height), cex=1.3,
          lab=expression(epsilon[u]*sigma*T[u]^4))
straightarrow(c(upper.emission.pos,upper.height-rady/2), 
              c(upper.emission.pos,lower.height+rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Radiation emitted from upper to space
straightarrow(c(upper.emission.pos,upper.height+rady/2), 
              c(upper.emission.pos,upper.height+2.*rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Radiation emitted from lower to surface
textplain(c(lower.emission.pos,lower.height), cex=1.3,
          lab=expression(epsilon[L]*sigma*T[L]^4))
straightarrow(c(lower.emission.pos,lower.height-rady/2), 
              c(lower.emission.pos,surface.height+rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Radiation from lower to upper
straightarrow(c(lower.emission.pos,lower.height+rady/2), 
              c(lower.emission.pos,upper.height-rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Radiation emitted from surface to lower
textplain(c(surface.emission.pos,surface.height), cex=1.3,
          lab=expression(sigma*T[s]^4))
straightarrow(c(surface.emission.pos,surface.height+rady/2), 
              c(surface.emission.pos,lower.height-rady), arr.pos=1,
              endhead=TRUE, lwd=arrow.thickness)

# Add temperature labels in white rectangles
textrect(c(temp.pos, upper.height), 0.10, rady=0.05, shadow.size=0, 
         lab=expression(T[L]== phantom(0) ), adj=c(1,0.5))
textrect(c(temp.pos, lower.height), 0.10, rady=0.05, shadow.size=0, 
         lab=expression(T[L]== phantom(0) ), adj=c(1.,0.5))
textrect(c(temp.pos, surface.height), 0.10, rady=0.05, shadow.size=0, 
         lab=expression(T[s]== phantom(0) ), adj=c(1.,0.5))