Overview

Simulation of a tidal sandflat environment with NicheMapR. The idea is to use the microclimate model with the shore option turned, and giving the model a vector of hourly rainfall such that it rains heavily when the tide is in.

The simulation is tested against Yoann Thomas’ data for a location in Senegal, West Africa. The data include observed sand (~2.5 cm) temperatures and water temperatures, as well as the presence/absence of the tide. 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(raster)
## Loading required package: sp
library(microclima)

loc <- c(-16.6798, 13.9236)

tzoff <- "+0"
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)

Sandflat observation data

The plot shows the sand temperature (red) and water temperature (blue) in relation to the tide presence (grey).

load("data_temp_Yoann.Rdata")
data$dates <- as.POSIXct(data$date, format = "%Y-%m-%d %H:%M:%S", 
    tz = tzone)

xlim <- c(as.POSIXct("2018-06-20"), as.POSIXct("2018-11-05"))
plot(data$dates + as.numeric(tzoff) * 3600, data$temp_water, 
    type = "l", col = "light blue", xlim = xlim, ylim = c(20, 
        40), xlab = "time", ylab = "temperature, deg C")
points(data$dates + as.numeric(tzoff) * 3600, data$water * 25, 
    col = "grey", type = "h")
points(data$dates + as.numeric(tzoff) * 3600, data$temp_sand, 
    type = "l", col = "red")

xlim <- c(as.POSIXct("2018-09-15"), as.POSIXct("2018-10-01"))
plot(data$dates + as.numeric(tzoff) * 3600, data$temp_water, 
    type = "l", col = "light blue", xlim = xlim, ylim = c(20, 
        40), xlab = "time", ylab = "temperature, deg C")
points(data$dates + as.numeric(tzoff) * 3600, data$water * 25, 
    col = "grey", type = "h")
points(data$dates + as.numeric(tzoff) * 3600, data$temp_sand, 
    type = "l", col = "red")

Run the microclimate model without a tidal effect

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 ‘sandflat data’.

Note that the depths to simlated are specified by DEP. You can change these depths, e.g. could make it 0, 2, 4, … depending on where you think the most representative logger depth is.

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

DEP <- c(0, 2.5, 5, 10, 15, 20, 30, 50, 100, 200)  # depths to simulate

# run the model note that 'save' is set to 2 because the
# model was run previously and the NCEP inputs saved in the
# folder 'sandflat data/'. If you haven't run it already,
# you'll need to set 'save = 1'
setwd("sandflat data/")
micro <- micro_ncep(DEP = DEP, loc = loc, dstart = dstart, dfinish = dfinish, 
    reanalysis = FALSE, save = 2, coastal = T, cap = 0, soilgrids = 1)
## loading SoilGrids data from previous run 
## 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 -16.6798 lat 13.9236 
##    user  system elapsed 
##   12.33    0.04   12.66
setwd("..")
metout <- as.data.frame(micro$metout)
soil <- as.data.frame(micro$soil)
metout <- cbind(micro$dates, metout)
soil <- cbind(micro$dates, soil)

xlim <- c(as.POSIXct("2018-06-20"), as.POSIXct("2018-11-05"))
plot(data$dates + as.numeric(tzoff) * 3600, data$temp_water, 
    type = "l", col = "light blue", xlim = xlim, ylim = c(20, 
        40), xlab = "time", ylab = "temperature, deg C")
points(data$dates + as.numeric(tzoff) * 3600, data$water * 25, 
    col = "grey", type = "h")
points(data$dates + as.numeric(tzoff) * 3600, data$temp_sand, 
    type = "l", col = "red")
points(micro$dates, soil$D2.5cm, type = "l")

xlim <- c(as.POSIXct("2018-09-15"), as.POSIXct("2018-10-01"))
plot(data$dates + as.numeric(tzoff) * 3600, data$temp_water, 
    type = "l", col = "light blue", xlim = xlim, ylim = c(20, 
        40), xlab = "time", ylab = "temperature, deg C")
points(data$dates + as.numeric(tzoff) * 3600, data$water * 25, 
    col = "grey", type = "h")
points(data$dates + as.numeric(tzoff) * 3600, data$temp_sand, 
    type = "l", col = "red")
points(micro$dates + as.numeric(tzoff) * 3600, soil$D2.5cm, type = "l")

Construct a tide table and hourly rainfall (tidal innundation) vector

This section builds the tides array, based on the ‘water’ column in the input data above and on the hourly water temperature data ‘temp_water’.

# aggregate tide presence/absence and water temperature by
# hour
agg_tide <- aggregate(data$water, by = list(format(data$dates, 
    "%Y-%m-%d %H")), FUN = max)
agg_twater <- aggregate(data$temp_water, by = list(format(data$dates, 
    "%Y-%m-%d %H")), FUN = mean)

# create tides vector
tides <- matrix(data = 0, nrow = length(dates), ncol = 3)
tides[, 1] <- 0  # default tide is out
tides[dates %in% as.POSIXct(agg_tide$Group.1, format = "%Y-%m-%d %H", 
    tz = tzone), 1] <- agg_tide$x  # tide in periods
tides[dates %in% as.POSIXct(agg_twater$Group.1, format = "%Y-%m-%d %H", 
    tz = tzone), 2] <- agg_twater$x  # water temperature

# set up a rain vector to simulate innundation by tide
microclima <- micro$microclima.out  # extract raw input from NCEP obtained via microclima
raindaily <- microclima$dailyprecip  # get daily rain
# construct an hourly rain vector
hourlyrain <- rep(0, length(dates))  # empty
hourlyrain[seq(1, length(hourlyrain), 24)] <- raindaily  # add actual rainfall at midnight
hourlyrain[dates %in% as.POSIXct(agg_tide$Group.1[agg_tide$x == 
    1], format = "%Y-%m-%d %H", tz = tzone)] <- 100  # add arbitrary rain when tide is in

Run the microclimate model simulating the tide effect

Now simulate the sandflat microclimate with the shore option turned on, the substrate set to sand as above, and simulated surface temperature (black), accounting for the effect of the tide coming in.

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 ‘sandflat data’.

Also note that the ERR parameter had to be set a bit higher and the snow model turned off because the model was having trouble converging on a solution.

The results are compared with 2.5 cm predicted sand temperatures.

rainhourly <- 1  # use hourly rainfall input
shore <- 1  # use shore model (i.e. tide table used to impose strong convection at water temperature when tide is in)
snowmodel <- 0  # turn this off - helps with stability of solution
ERR <- 2.5  # integrator error - has to be set higher than default of 1.5
maxpool <- 1  # maximum depth of water pooling on surface - make small because tide simulated by rainfall
setwd("sandflat data/")
micro <- micro_ncep(DEP = DEP, rainhourly = rainhourly, rainhour = hourlyrain, 
    shore = shore, snowmodel = snowmodel, ERR = ERR, tides = tides, 
    loc = loc, dstart = dstart, dfinish = dfinish, reanalysis = FALSE, 
    save = 2, coastal = T, cap = 0, soilgrids = 1, maxpool = 1)
## loading SoilGrids data from previous run 
## 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 -16.6798 lat 13.9236 
##    user  system elapsed 
##   35.32    0.03   37.61
setwd("..")
metout <- as.data.frame(micro$metout)
soil <- as.data.frame(micro$soil)
metout <- cbind(micro$dates, metout)
soil <- cbind(micro$dates, soil)

# plot the final results

# function for plotting (note predicted value at 2.5 cm is
# chosen to plot)
plot_data <- function(xlim) {
    plot(data$dates + as.numeric(tzoff) * 3600, data$temp_water, 
        type = "l", col = "light blue", xlim = xlim, ylim = c(20, 
            40), xlab = "time", ylab = "temperature, deg C")
    points(data$dates + as.numeric(tzoff) * 3600, data$water * 
        25, col = "grey", type = "h")
    points(data$dates + as.numeric(tzoff) * 3600, data$temp_water, 
        type = "l", col = "light blue")
    points(data$dates + as.numeric(tzoff) * 3600, data$temp_sand, 
        type = "l", col = "red")
    points(micro$dates, soil$D2.5cm, type = "l")
}

# whole time series
xlim <- c(as.POSIXct("2018-06-20"), as.POSIXct("2018-11-05"))
plot_data(xlim)

# two-weekly plots
xlim <- c(as.POSIXct("2018-06-15"), as.POSIXct("2018-07-01"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-07-01"), as.POSIXct("2018-07-15"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-07-15"), as.POSIXct("2018-08-01"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-08-01"), as.POSIXct("2018-08-15"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-08-15"), as.POSIXct("2018-09-01"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-09-01"), as.POSIXct("2018-09-15"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-09-15"), as.POSIXct("2018-10-01"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-10-01"), as.POSIXct("2018-10-15"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-10-15"), as.POSIXct("2018-11-01"))
plot_data(xlim)

xlim <- c(as.POSIXct("2018-11-01"), as.POSIXct("2018-11-15"))
plot_data(xlim)

An aside - how soil properties are specified in NicheMapR

The code below shows how you’d set substrate hydraulic properties manually as well as how the ‘SoilGrids’ database https://soilgrids.org/#!/?layer=ORCDRC_M_sl2_250m&vector=1 is used within micro_ncep.

# here's how you'd set the set substrate properties, but for
# now using the SoilGrids database
BD <- rep(1.3, 19)  # bulk density
PE <- BD  # dummy vector
BB <- BD  # dummy vector
KS <- BD  # dummy vector
soiltype <- 1  # sand
PE[1:19] <- CampNormTbl9_1[soiltype, 4]  # fill air entry potential vector
BB[1:19] <- CampNormTbl9_1[soiltype, 5]  # fill Campbell's B parameter vector
KS[1:19] <- CampNormTbl9_1[soiltype, 6]  # fill saturated hydraulic conductivity vector

# code internal to micro_ncep used to get soil properties
# from soilgrids
DEP <- c(0, 2.5, 5, 10, 15, 20, 30, 50, 100, 200)
x <- t(as.matrix(as.numeric(c(loc[1], loc[2]))))
ov <- jsonlite::fromJSON(paste0("https://rest.soilgrids.org/query?lon=", 
    x[1], "&lat=", x[2], ",&attributes=BLDFIE,SLTPPT,SNDPPT,CLYPPT"), 
    flatten = TRUE)
soilpro <- cbind(c(0, 5, 15, 30, 60, 100, 200), unlist(ov$properties$BLDFIE$M)/1000, 
    unlist(ov$properties$SLTPPT$M), unlist(ov$properties$SNDPPT$M), 
    unlist(ov$properties$CLYPPT$M))
colnames(soilpro) <- c("depth", "blkdens", "clay", "silt", "sand")
# Now get hydraulic properties for this soil using Cosby et
# al. 1984 pedotransfer functions.
soil.hydro <- pedotransfer(soilpro = as.data.frame(soilpro), 
    DEP = DEP)
PE <- soil.hydro$PE
BB <- soil.hydro$BB
BD <- soil.hydro$BD
KS <- soil.hydro$KS