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.

Consider a simple planet with two isothermal layers of graybody atmosphere in radiative equilibrium with a blackbody surface.

Assume the atmosphere is transparent to solar radiation, the lower atmospheric layer has broadband thermal emissivity denoted by \(\epsilon_1\), and the upper atmospheric layer has broadband thermal emissivity denoted by \(\epsilon_2\). Let the planetary albedo be \(\alpha\), and the solar flux at teh top of the atmosphere be designated \(S_0\).

The temperature of the surface, the lower atmospheric layer, and the upper atmospheric layer are denoted \(T_S\), \(T_1\), and \(T_2\) respectively.

Assume that all nonradiative energy fluxes (e.g., turbulence, convection, evaporation, condensation) can be represented as generalized net “convective” heat fluxes \(H_1\) (from the surface to the lower atmosphere) and \(H_2\) (from the lower to the upper atmosphere. (Of course, there can’t be any convective fluxes out the top of the upper atmosphere!)

All emission is isotropic. By the Stefan-Boltzmann Law, the surface emits thermal infrared radiation at the rate \(\sigma T_S^4\). The lower atmospheric layer emits energy both upward and downward at the rate \(\epsilon_1 \sigma T_1^4\), and the upper atmospheric layer emits energy at the rate \(\epsilon_2 \sigma T_2^4\).

By Kirchoff’s Law the absorptivity of each atmospheric layer in thermal infrared wavelengths is precisely the same as its emissivity. Since the surface is a blackbody, it absorbs all radiation incident on it.

The surface absorbs energy from the Sun at a rate \(\frac{S_0}{4}(1-\alpha)\). It also absorbs energy at a rate \(\epsilon_1 \sigma T_1^4\) from the lower atmospheric layer. Since a fraction \(\epsilon_2\) of the longwave radiation emitted downward by the upper atmospheric layer is absorbed by the lower atmospheric layer, the surface absorbs this radiation at a rate \((1-\epsilon_1)\epsilon_2 \sigma T_2^4\).

The lower atmospheric layer is transparent to solar radiation, but absorbs a fraction \(epsilon_1\) of the longwave radiation incident upon it. From the surface, it absorbs upwelling energy at the rate \(\epsilon_1 \sigma T_S^4\), and from the upper layer it absorbs downwelling energy at the rate \(\epsilon_1 \epsilon_2 \sigma T_2^4\).

The upper atmospheric layer is also transparent to solar radiation, but absorbs a fraction \(\epsilon_2\) of the longwave radiation incident upon it. Since a fraction \(\epsilon_1\) of the longwave radiation emitted upward by the surface is absorbed by the lower atmospheric layer (see previous paragraph), the upper layer absorbs upwelling energy emitted from the surface at the rate \(\epsilon_2 (1-\epsilon_1) \sigma T_S^4\). From the lower atmospheric layer, the upper layer absorbs upwelling energy at the rate \(\epsilon_1 \epsilon_2 \sigma T_1^4\).

Using the assumption that the surface and each atmospheric layer is in radiative equilibrium, we can now write a radiation flux budget for the surface, for each atmospheric layer, and for the planet as a whole (defined at the top of the atmosphere). In each of the energy balances that follow, sources of energy for the layer are grouped on the left-hand sides, and losses to the right of the equal signs:

\[ \begin{array}{lrclr} Surface: & \frac{S_0}{4}(1-\alpha) + \epsilon_1 \sigma T_1^4 + (1-\epsilon_1) \epsilon2 \sigma T_2^4 &=& \sigma T_S^4 + H_S & (1)\\ Lower Layer: & \epsilon_1 \sigma T_S^4 + \epsilon_1 \epsilon_2 \sigma T_2^4 + H_S &=& 2 \epsilon_1 \sigma T_1^4 + H_L & (2)\\ Upper Layer:& (1-\epsilon_1) \epsilon_2 \sigma T_S^4 + \epsilon_1 \epsilon_2 \sigma T_1^4 + H_L &=& 2 \epsilon_2 \sigma T_2^4 & (3)\\ \end{array} \]

The energy budgets (1) through (3) constitute a system of three independent equations in three unknowns: \(T_S\), \(T_1\), and \(T_2\).

Adding equations (1) through (3) provides a statement that for the planet as a whole the absorbed solar flux is balanced by the total outgoing longwave radiation (OLR) at the top of the atmosphere (TOA):

\[ \begin{array}{lrclr} Planet (TOA):& \frac{S_0}{4}(1-\alpha) &=& (1-\epsilon_1) (1-\epsilon_2) \sigma T_S^4 + \epsilon_1 (1-\epsilon_2) \sigma T_1^4 + \epsilon_2 \sigma T_2^4. & (4)\\ \end{array} \]

We can rearrange equations (1) through (3) to obtain an equivalent matrix equation:

\[ \begin{bmatrix} -1 & \epsilon_1 & (1-\epsilon_1)\epsilon_2 \\ \epsilon_1 & -2 \epsilon_1 & \epsilon_1 \epsilon_2 \\ (1-\epsilon_1) \epsilon_2 & \epsilon_1 \epsilon_2 & -2 \epsilon_2 \end{bmatrix} \times \left[ \begin{array}{c} T_S^4 \\ T_1^4 \\ T_2^4 \end{array} \right] = \frac{1}{\sigma} \left[ \begin{array}{c} H_s-\frac{S_0}{4}(1-\alpha) \\ H_L-H_s \\ -H_L \end{array} \right]. \]

To solve, we multiply the inverse of the square matrix by the vector on the right-hand side and then take the fourth root of each element of the resulting vector to obtain each temperature.

This is not as complicated as it sounds. To see how easy this is to set up and solve in R, see the “Website Code” tab to the right!

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!

**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'))
)
)
))
```

**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)
})
```

**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))
}
```

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))
```