Overview

Simulation of a tidepool thermal environment with NicheMapR. The idea is to use the microclimate model with the shore option turned, and then the ectotherm model with the container option turned on.

First the microclimate model is run with a table of tidal conditions, tides, is provided. This table has presence/absence (1 or 0) in the first column indicating if the tide is in at the location, and the water temperature in the second column (the third column indicates if wave spash, i.e. surface evaporation, is present).

Second, the ectotherm model is used to simulate a container sunk into the ground of a specified depth and diameter using the transient heat budget calculations of the ectotherm model to work out mean water temperature. The container is simulated to be full when the tide is in but will start drying out and changing temperature as soon as the tide is out, as a function of the microclimatic conditions.

The simulation is tested against Brian Helmuth et al.’s tidepool data for a location in Maine USA. The data include observed tidepool temperatures and water temperatures, as well as the tidal height of the pools. The forcing data are obtained from NCEP via the micro_ncep function that integrates with the microclima package.

Setup and location

Load the libraries and get the digital elevation model for the site via package elevatr and microclima.

library(NicheMapR)
library(microclima)
library(raster)
## Loading required package: sp
# Shapes are estimates but we can probably get more exact
# measurements Large tide pool: max depth 1.07m; shape
# 'round'; 2.5m greatest length small tide pool: max depth
# 0.22m; shape 'oval'; 0.83m greatest length Latitude:
# 42.427494 Longitude: -70.920384

tzoff <- -0  # need to check carefully how time zones work out given GMT for micro_ncep output
loc <- c(-70.920384, 42.427494)  # Maine
tzone <- paste("Etc/GMT", tzoff, sep = "")

# plot DEM
dem <- get_dem(lat = loc[2], long = loc[1])
## Merging DEMs
## Reprojecting DEM to original projection
## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
##  +init=epsg:3395 +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
plot(dem)
e <- extent(dem)
xy <- data.frame(x = (e@xmin + e@xmax)/2, y = (e@ymin + e@ymax)/2)
coordinates(xy) = ~x + y
points(xy, pch = 16)

Tidepool data input

# read in tidepool data
tp_large <- read.csv("04132018_tidepool_largewithtide.csv", stringsAsFactors = FALSE, 
    skip = 4, header = FALSE)
colnames(tp_large) <- c("#", "day", "time", "date", "low", "high", 
    "temperature", "tideheight")
tp_large$dates <- as.POSIXct(paste(tp_large$day, tp_large$time), 
    format = "%m/%d/%y %H:%M", tz = tzone)

tp_small <- read.csv("04132018_tidepool_smallwithtide.csv", stringsAsFactors = FALSE, 
    skip = 4, header = FALSE)
colnames(tp_small) <- c("#", "day", "time", "date", "low", "high", 
    "temperature", "tideheight")
tp_small$dates <- as.POSIXct(paste(tp_small$day, tp_small$time), 
    format = "%m/%d/%y %H:%M", tz = tzone)

# read in water temperature
twater <- read.csv("DailyWaterTemp.csv", stringsAsFactors = FALSE)
colnames(twater) <- c("date", "twater", "stdev")
twater$dates <- as.Date(twater$date, "%m/%d/%y", tz = tzone)

# define tidal periods
tp_large$tide <- 0
tp_large$tide[tp_large$tideheight > 0.7] <- 1  # 0.7 m tidal elevation
# aggregate to hourly
tp_small_agg_tide <- aggregate(tp_small$tide, by = list(format(tp_small$dates, 
    "%Y-%m-%d %H")), FUN = max)

tp_small$tide <- 0
tp_small$tide[tp_small$tideheight > 1.1] <- 1  # 1.1 m tidal elevation
# aggregate to hourly
tp_small_agg_tide <- aggregate(tp_small$tide, by = list(format(tp_small$dates, 
    "%Y-%m-%d %H")), FUN = max)

par(mfrow = c(2, 2))
# plot small tidepool data
plot(tp_large$dates, tp_large$temperature, type = "l")
plot(tp_large$dates, tp_large$low, type = "l", col = "red")
points(tp_large$dates, tp_large$high, type = "l", col = "blue")

# plot small tidepool data
plot(tp_small$dates, tp_small$temperature, type = "l")
plot(tp_small$dates, tp_small$low, type = "l", col = "red")
points(tp_small$dates, tp_small$high, type = "l", col = "blue")

Construct a tide table

This section builds the tides array, based on the tide in/out data above and on an hourly spline of the daily water temperature data. The plot shows the water temperature (red) in relation to the tide presence (grey).

# set simulation period
dstart <- "01/02/2018"
dfinish <- "30/11/2018"
dates <- seq(ISOdate(2018, 2, 1) - 3600 * 12, ISOdate(2018, 12, 
    1) - 3600 * 13, by = "hours", tz = "UTZ")

# create a tide table for NicheMapR
tides <- matrix(data = 0, nrow = length(dates), ncol = 3)
tides[, 1] <- 0  # default is tide out
# now sub in known periods when tide was in on the basis of
# date
tides[dates %in% as.POSIXct(tp_small_agg_tide$Group.1, format = "%Y-%m-%d %H", 
    tz = tzone), 1] <- tp_small_agg_tide$x

# spline daily water temperatures
rng <- range(as.POSIXct(twater$dates))  # date range
tt <- seq(rng[1], rng[2], by = "hour")  # time sequence to spline to
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
z <- spline(cbind(as.POSIXct(twater$dates), twater$twater), xout = tt)  # spline
twater2 <- cbind(as.data.frame(tt), z$y)  # add dates
par(mfrow = c(1, 1))
plot(twater2, type = "l")

colnames(twater2) <- c("date", "twater")

# add the splined daily water temps to the tide table
tides[which(dates == min(twater2$date - tzoff * 3600)):which(dates == 
    max(twater2$date - tzoff * 3600)), 2] <- twater2$twater

# plot the tide data thus far

xlim <- c(as.POSIXct("2018-04-15"), as.POSIXct("2018-07-05"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(dates, tides[, 1] * 20, col = "grey", type = "h")
points(dates, tides[, 2], col = "blue", type = "l")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-05-20"), as.POSIXct("2018-06-05"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(dates, tides[, 1] * 20, col = "grey", type = "h")
points(dates, tides[, 2], col = "blue", type = "l")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

Run the microclimate model simulating the tide effect

Now simulate the rock platform microclimate with the shore option turned on, the substrate set to rock (Density is the same as BulkDensity and cap off). The plot now includes solar radiation (orange) and simulated surface temperature (black), accounting for the effect of the tide coming in on the the rock platform.

Note that this microclimate model run assumes it has been run previously with save = 1 and that the forcing data are saved in a subdirectory called ‘tidepool data’.

# run microclimate model with tide on, for rock

# note that 'save' is set to 2 because the model was run
# previously and the NCEP inputs saved in the folder
# 'tidepool data/'. If you haven't run it already, you'll
# need to set 'save = 1'
setwd("tidepool data/")
micro <- micro_ncep(tides = tides, loc = loc, dstart = dstart, 
    runmoist = 0, snowmodel = 0, dfinish = dfinish, reanalysis = FALSE, 
    write_input = 0, writecsv = 0, save = 2, BulkDensity = 2.56, 
    shore = 1, coastal = T, cap = 0)
## loading met data from previous run 
## running microclimate model for  303  days from  2018-02-01  to  2018-11-30 23:00:00  at site  long -70.920384 lat 42.427494 
##    user  system elapsed 
##    9.70    0.03   10.80
setwd("..")

# extract results and add dates
metout <- as.data.frame(micro$metout)
soil <- as.data.frame(micro$soil)
metout <- cbind(micro$dates, metout)
soil <- cbind(micro$dates, soil)

# plot again
xlim <- c(as.POSIXct("2018-04-15"), as.POSIXct("2018-07-05"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, soil$D0cm, col = "black", type = "l")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-05-20"), as.POSIXct("2018-06-05"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, soil$D0cm, col = "black", type = "l")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

Simulate tidepool with ectotherm container model

Finally, the ectotherm model is run with container set to 1 to turn on the container model. A vector of daily ‘rainfall’ values is constructed so that it is raining arbitrarily heavily whenever the tide is in. In addition, make air and sky temperature equal water temperature when the tide is in and make the wind speed arbitrarily large when tide comes in. Note that the calculations with the ectotherm model running like this take a while (approx 2 mins on my machine to run this time series of just under a year).

# create a vector of 'rainfall'
poolfill <- aggregate(tides[, 1], by = list(format(micro$dates, 
    "%Y-%m-%d")), FUN = max)$x * 50000
poolfillhr <- tides[, 1] * 50000

# set tide microclimate make air temp water temp when tide is
# in
micro$metout[which(tides[, 1] == 1), 3] <- tides[which(tides[, 
    1] == 1), 2]
# make sky temp water temp when tide is in
micro$metout[which(tides[, 1] == 1), 14] <- tides[which(tides[, 
    1] == 1), 2]
# make windspeed high when tide is in
micro$metout[which(tides[, 1] == 1), 7] <- micro$metout[which(tides[, 
    1] == 1), 7] * 5000
# make RH 100% when tide is in
micro$metout[which(tides[, 1] == 1), 5] <- 100
# make solar 0 when tide is in
micro$metout[which(tides[, 1] == 1), 13] <- 0


# run the simulation
pooldepth <- 200  # pool maximum depth, mm
poolwidth <- 800  # pool width (diameter), mm
alpha_max <- 0.9  # pool solar absorptivity
alpha_min <- alpha_max
ecto <- ectotherm(alpha_max = alpha_max, alpha_min = alpha_min, 
    container = 1, conth = pooldepth, contw = poolwidth, live = 0, 
    rainfall = poolfill, rainhr = poolfillhr, contwet = 100, 
    contonly = 1, continit = 0, contype = 1)
## running ectotherm model ...
## runtime 144.97 seconds
environ <- as.data.frame(ecto$environ)
enbal <- as.data.frame(ecto$enbal)
masbal <- as.data.frame(ecto$masbal)
debout <- as.data.frame(ecto$debout)

# plot simulated pool depth
xlim <- c(as.POSIXct("2018-04-15"), as.POSIXct("2018-07-05"))
plot(micro$dates, environ$CONDEP, type = "l", xlab = "date", 
    ylab = "pool depth, mm")

plot(micro$dates, environ$WATERTEMP, type = "l")

# plot again
xlim <- c(as.POSIXct("2018-04-15"), as.POSIXct("2018-07-05"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$WATERTEMP, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-04-15"), as.POSIXct("2018-05-01"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$TC, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-05-01"), as.POSIXct("2018-05-15"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$TC, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-05-15"), as.POSIXct("2018-06-01"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$TC, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-06-01"), as.POSIXct("2018-06-15"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$TC, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-06-15"), as.POSIXct("2018-07-01"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$TC, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")

xlim <- c(as.POSIXct("2018-07-01"), as.POSIXct("2018-07-15"))
plot(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 1] * 20, col = "grey", type = "h")
points(micro$dates, metout$SOLR/100, col = "orange", type = "h", 
    lwd = 1.5)
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")
points(micro$dates, tides[, 2], col = "blue", type = "l")
points(micro$dates, environ$TC, type = "l", col = "black")
points(tp_small$dates, tp_small$temperature, type = "l", ylim = c(0, 
    30), xlim = xlim, col = "red")