Cascading Pools with Steady Unit Inflow

Set drainage time for each pool

Pools and Flows Over Time

MODEL DESCRIPTION

Demonstration of flows through a cascading series of pools with different “drain rates.” At the top of the cascade is an inflow (rate = 1). Each pool then drains at a rate that’s directly proportional to the volume of the pool.

For each pool, \[\frac{dP}{dt} = inflow - \frac{P}{\tau}\] where P is the volume in the pool and \(\tau\) is the e-folding timescale for the pool to drain downstream.

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 (run.pools.R) that integrates the cascading pools model
  • A script (plot.pools.R) that plots the output
  • 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!

run.pools.R:

run.pools <- function(time1=3, time2=10, time3=100){

  # Create drainage and flow transfer matrix
  drain <- c(1./time1, 1./time2, 1./time3)
  flow <- matrix(c(-drain[1],     0,        0,
                    drain[1], -drain[2],    0,
                      0,       drain[2], -drain[3]),
                   byrow=TRUE, nrow=3)

  # Inflow into top pool only
  inflow <- c(1,0,0)

  # Compute steady-state solution
  steady <- solve(flow, -inflow)

  # Initialize pools & outflows to zeros 
  nTime <- 2*max(c(time1, time2, time3))
  pool <- matrix(replicate(3*nTime, 0.), nrow=nTime)
  outflow <- pool

  # Run the model forward in time
  for (i in 1:(nTime-1)) {
    pool[i+1, ] <- pool[i, ] + inflow + flow %*% pool[i, ]
    outflow[i+1,] <- drain * pool[i+1, ]
  }

  # return timeseries of each pool to calling program
  return(list(pool=pool, outflow=outflow, steady=steady))  

}

plot.pools.R:

plot.pools <- function(pools){

  # Set up plot parameters for a 2-panel (landscape) page
  orig.par <- par(no.readonly=TRUE)  # Remember changable parameters to reset later
  par(mfrow=c(3,2), cex=1.2, mar=c(2,2,2,2))  

  # Top row shows volume and outflow from top pool
  plot(pools$pool[,1], type='l', col='red', lwd=4, main='Top Volume',
       ylim=c(0,pools$steady[1]), ylab='', xlab='time')
  abline(h=pools$steady[1], col='blue', lty=2, lwd=4)
  plot(pools$outflow[,1], type='l', col='blue', lwd=4, main='Top outflow',
       ylim=c(0,1), ylab='', xlab='time') 

  # Middle row shows volume and outflow from middle pool
  plot(pools$pool[,2], type='l', col='red', lwd=4, main='Middle Volume',
       ylim=c(0,pools$steady[2]), ylab='', xlab='time')
  abline(h=pools$steady[2], col='blue', lty=2, lwd=4)
  plot(pools$outflow[,2], type='l', col='blue', lwd=4, main='Middle outflow',
       ylim=c(0,1), ylab='', xlab='time') 

  # Bottom row shows volume and outflow from bottom pool
  plot(pools$pool[,3], type='l', col='red', lwd=4, main='Bottom Volume',
       ylim=c(0,pools$steady[3]), ylab='', xlab='time')
  abline(h=pools$steady[3], col='blue', lty=2, lwd=4)  
  plot(pools$outflow[,3], type='l', col='blue', lwd=4, main='Bottom outflow',
       ylim=c(0,1), ylab='', xlab='time') 

  par(orig.par)

}

This is the user interface that actually controls this website

ui.R:

library(shiny)
library(markdown)
library(knitr)

# Define user interface for logistic growth model
shinyUI(pageWithSidebar(

  #  Application title
  headerPanel("Logistic Growth with Nutrient Limitation"),

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

    img(src='growth.jpg'),

    h4('Set Models Parameter Below'),            

    # Set drainage times for each pool
    sliderInput("max.growth", "Maximum Growth Rate (percent)", 
                min = 1, max = 100, value = 10, step= 1),
    sliderInput("lifetime", "Average Lifetime (time steps)", 
                min = 1, max = 50, value = 10, step= 1),
    sliderInput("limit", "Nutrient-Limited Maximum Population", 
                min = 1, max = 1000, value = 100, step=1),
    sliderInput("nTime", "Number of time steps", 
                min = 10, max = 500, value = 100, step=10)
  ),

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

      # First tab shows plots of model output for both control and modified atmosphere
      tabPanel("Output", 
               h4('Population Over Time'), 
               p('  death <- max.growth / lifetime'),
               p('  growth[i+1] <- max.growth * (1 - life[i] / limit)'),
               p('  life[i+1] <- life[i] * (1 + growth[i] - death)'),
               plotOutput('population.plot')), 

      # Second tab displays a brief model description
      tabPanel("Model Description", 
               includeHTML('doc/model.description.html')),

      # Third tab displays website code
      tabPanel("Website code", 
               includeMarkdown('doc/website.code.md'))
    )
  )
))

server.R:

library(shiny)

# Source required R scripts
source('logistic.R')        

shinyServer(function(input, output) {

  #  plot of model output
  output$population.plot <- renderPlot(
      plot(logistic(max.growth=input$max.growth/100, lifetime=input$lifetime, 
                    limit=input$limit, nTime=input$nTime),
           type='l', lwd=10, col='darkgreen', ylab='', xlab='time', 
           cex.lab=2, cex.axis=2)
  )
})