Sunshine on a Perfectly Cloudless Day

Where and When

Diurnal Cycle of Insolation

MODEL DESCRIPTION

This is a very simple model of incoming solar radiation at the Earth’s surface, intended to be used for teaching and learning about weather and climate.

The user specifies the season (month and day), the latitude, and the pressure. The program then uses simple spherical geometry to calculate the incident solar radiation on a horizontal plane, accounting for atmospheric attentuaion by absorption and scattering using a very simple expoential function suitable for clear skies. There’s no way to account for clouds in this simple model!

The model follows the discussion in Bonan (Chapter 4). in the absence of atmospheric scattering and absorption, the power of the Sun’s radiation incident on a horizontal surface is given by Bonan’s equation 4.4:

\[\begin{aligned} S_H = \frac{S_C}{r_v^2}cos Z \end{aligned}\]

where \(S_C\) = 1367 W m-2 is the “solar constant,” \(r_v\) is the “radius vector” (the ratio of the actual distance between the Sun and the Earth to its annual average), and \(cosZ\) is the cosine of the solar zenith angle. When the cosine of the solar zenith angle is negative, the Sun is below the horizon and we can just set the incoming solar radiation to zero.

The cosine of the solar zenith angle is given by Bonan’s equation 4.1:

\[ \begin{aligned} cos Z = sin \phi sin \delta + cos \delta cos \phi cos h \end{aligned} \]

where is latitude, \(\delta\) is declination (the latitude where the Sun is directly overhead at noon), and \(h\) is the hour angle (\(h = 0\) at noon, \(h = \pi/2\) at 6 PM, \(h = \pi\) at midnight, and \(h = 3\pi/2\) at 6 AM).

To calculate the hour angle, simply compute

\[ \begin{aligned} h = -\pi + \frac{LST}{24} 2 \pi \end{aligned} \]

where LST is local solar time (or local standard time, same thing).

Declination (latitude of the sub-Solar point) is given by

\[ \begin{aligned} \delta = \epsilon sin \theta \end{aligned} \]

where \(\epsilon\) is the tilt of the Earth’s axis (currently \(\epsilon = 0.4091\) radians) and \(\theta\) is the angular (seasonal) position of Earth in it’s orbit. The seasonal position angle \(\theta\) is measured from the vernal equinox (March 21), so that \(\theta = \pi/2\) at the June solstice, \(\theta = \pi\) at the autumnal equinox (Sept 21), and \(\theta = 3\pi/2\) at the December solstice.

The radius vector (normalized Earth-Sun distance) depends on both the position of the Earth in its orbit and the “shape” of the orbit itself according to:

\[ \begin{aligned} \frac{1}{r_v} = \frac{R_0}{R_E} = \frac{(1-e^2)}{[1+e cos(\theta - \omega)]} \end{aligned} \]

where \(R_0\) is the average Sun-Earth distance, \(R_E\) is the actual Sun-Earth distance, \(e\) is the eccentricity of Earth’s orbit (deviation from a perfect circle), \(\theta\) is the seasonal position angle (see above), and \(\omega\) is the angle between the orbital position of the Earth’s perihelion and the vernal equinox. The “Milankovitch” orbital parameters \(\epsilon\), \(e\), and \(\omega\) vary quasi-periodically over geologic time due to the gravitational influences of the other planets. Present-day values are \(\epsilon\) = 0.4091, \(e\) = 0.0167 and \(\omega\) = 1.7963. See http://en.wikipedia.org/wiki/Insolation for a derivation and discussion.

The solar beam is also attenuated by atmospheric scattering and absorption. In a clear sky, incoming solar radiation at the Earth’s surface can be approximated using Bonan’s equation 4.6 as:

\[ \begin{aligned} S = S_H \tau^m \end{aligned} \]

where \(S_H\) is the incident solar radiation on a horizontal surface above the atmosphere (defined above), \(\tau\) is the atmospheric transmissivity, and \(m\) is the “air mass”" through which the solar beam must pass. For this problem, assume \(\tau\) = 0.7. The air mass is given by Bonan’s eqn 4.7:

\[ \begin{aligned} m = \frac{1}{cos Z} \frac{P}{P_S} \end{aligned} \]

where \(cos Z\) is the cosine of the solar zenith angle (defined above), \(P\) is atmospheric pressure and \(P_S\) = 1013.25 mb is the atmospheric pressure at sea level.

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:

  • A short program that calculates the radiation (insolation.R)
  • An even shorter program that plots the output (plot.insolation.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 module that calculates insolation for a given day and latitude

insolation.R:

insolation <- function(lat=40., month=6, day=22, press=1013.25){

  #######################
  # Given a location and date, give diurnal cycle ofincoming solar radiation 
  # at surface for each hour of the day. Return an array of 24 hourly values

  # For details and a derivation see Bonan chapter 4
  # Also some good info at http://en.wikipedia.org/wiki/Insolation
  #######################

    # Model parameters
    S0=1367         # Solar constant (Watts/m2)
    ecc=0.0167      # Eccentricity of Eart's orbit
    eps=0.4091      # Obliquity (tilt) of Eart's axis (radians)
    perih=1.7963    # Angle between March equinox & perihelion (radians)
    equinox <- 81   # Day of the year when March equinox occurs
    slp <- 1013.25  # sea-level pressure (mb)

    # find the day of the year
    # nDays is the number of days in each month
    nDays <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) 
    daynum <- 0  # count days starting with zero
    for (i in 1:(month-1)) daynum <- daynum + nDays[i]
    daynum <- daynum + day

    # Theta is the position angle of the Earth in its orbit, 
    # counted in radians starting at the Vernal Equinox (March 22)
    #     theta=0 is the March Equinox, theta=90 is June Solstice
    #     theta=180 is September Equinox, theta=270 is December Solstice)
    theta <- 2 * pi * (daynum - equinox) / 365

    # Convert latitude from degrees to radians
    lat <- lat * pi / 180.

    # Declination for this time of year (theta)
    # declination is the latitude where the Sun is directly overhead @ noon
    dec <- eps * sin(theta)

    # Earth-Sun distance (ratio of mean to actual)
    R0.div.RE <- (1 - ecc^2) / (1 + ecc * cos(theta - perih))

    # Create output array for hourly insolation values
    insol <- rep(NA,24)

    # Loop through the hours of the day
    for (LST in 0:23) {

      # Calculate insolation at surface
      h <- -pi + 2 * pi * LST/24  # Hour angle in radians of Sun past noon
      cos.z <- sin(lat)*sin(dec) + cos(lat)*cos(dec)*cos(h) # Bonan eq 4.1

      # Top of atmosphere insolation:  Bonan eq 4.4
      # (also see http://en.wikipedia.org/wiki/Insolation)
      toa <- max(S0 * R0.div.RE^2 * cos.z, 0)  # zero if sun is below horizon!

      # Attenuation of solar radiation by the atmosphere
      air.mass <- 1/cos.z * press/slp           # Bonan eq 4.7
      insol[LST+1] <- 1.1 * toa * 0.7^air.mass  # insol eq 4.6 +10% diffuse

  }

  return(insol)
}

This is module that plots the output of the radiation calculations

plot.insolation.R:

plot.insolation <- function(insol, lat=lat, month=month, day=day, press=press){

  # Make a plot of the diurnal cycle of incoming solar raditaion at the surface

  # Convert month number to month name
  month.name <- c('January','February','March','April','May','June',
                  'July','August','September','October','November','December')

  # Set plot margins
  bottomMargin <- 4
  leftMargin <- 5
  topMargin <- 1
  rightMargin <- .5

  # Set up plot margins
  orig.par <- par(no.readonly=TRUE)  # Remember changable parameters to reset later
  par(mar=c(bottomMargin,leftMargin,topMargin,rightMargin), oma=c(0,2.,1.,0))

  # Plot the diurnal cycle of radiation
  plot(0:23, insol, typ='l', col='red', lwd=5, xaxt='no',
       main = paste('Latitude =',lat,'on',month.name[month],day),
       xlab = 'Local Standard Time', ylab = expression(Watts~~m^{-2}),
       cex.main=1.4, cex.lab=1.3, cex.axis=1.2)
  axis(1, at=seq(0,24,3), cex.axis=1.2)

  # Add text for pressure and maximum insolation value
  text(3, 0.8*max(insol), paste('pressure =',press,'mb'))
  text(3, 0.7*max(insol), paste('max insolation =',format(max(insol), digits=4)))

  # 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 calculation of incoming solar radiation at surface
shinyUI(pageWithSidebar(

  #  Application title
  headerPanel("Sunshine on a Perfectly Cloudless Day"),

  # Sidebar with sliders that control various available options
  sidebarPanel(
    img(src='model.png'),
    h4('Where and When'),            

    selectInput("month", "Month:",
                c("January" = 1, 
                     "February" = 2,
                     "March" = 3,
                     "April" = 4,
                     "May" = 5,
                     "June" = 6,
                     "July" = 7,
                     "August" = 8,
                     "September" = 9,
                     "October" = 10,
                     "November" = 11,
                     "December" = 12),
    selected=6, selectize=FALSE),

    sliderInput("day", "Day:", 
                min = 1, max = 31, value = 15, step=1),

    sliderInput("latitude", "Latitude:", 
                min = -90, max = 90, value = 40, step=1),

    sliderInput("press", "Pressure (mb):",
                min=300, max=1020, value=1015, step=10)
  ),

   # Main panel consists of a tabset displays model output or model description
   mainPanel(
     tabsetPanel(
       tabPanel("Insolation", 
                h3("Diurnal Cycle of Insolation"),
                plotOutput("insolation.plot"),
                tags$style(type="text/css", ".tab-content { overflow: visible; }")), 
       tabPanel("Description", 
                includeHTML('doc/model.description.html')),
       tabPanel("Website Code",
                includeMarkdown('doc/website.code.md'))
     )
  )
))  

Here is the server module that responds to user interface events

server.R:

library(shiny)

# Source required R scripts
source('model/insolation.R')        
source('model/plot.insolation.R')        

shinyServer(function(input, output) {

  # Table of all the values in modelOutput
  output$insolation.plot <- renderPlot(
    plot.insolation(insolation(lat=input$latitude, 
                               month=as.numeric(input$month), 
                               day=input$day,
                               press=input$press), 
                  lat=input$latitude, month=as.numeric(input$month), 
                  day=input$day, 
                  press=input$press)
  )
}
)