micro_era5: ERA5 implementation of the microclimate model, assuming ERA5...

In mrke/NicheMapR: R implementation of Niche Mapper software for biophysical modelling

ERA5 implementation of the microclimate model, assuming ERA5 grids have been downloaded with package mcera5, and using package microclima downscaling for local topography.

Description

An implementation of the NicheMapR microclimate model that integrates the ERA5 hourly weather data and the elevatr package for obtaining DEM using downscaling functions from the microclima package, largely following the methods described in Kearney, M. R., Gillingham, P. K., Bramer, I., Duffy, J. P., & Maclean, I. M. D. (2019). A method for computing hourly, historical, terrain-corrected microclimate anywhere on Earth. Methods in Ecology and Evolution.

Usage

micro_era5(loc = c(-91.415669, -0.287145), dstart = "01/01/2019", dfinish = "31/07/2019",
REFL = 0.15, slope = 0, aspect = 0, DEP = c(0, 2.5,  5,  10,  15,  20,  30,  50,  100,  200), minshade = 0, maxshade = 90,
Usrhyt = 0.01, ...)

Arguments

loc	Longitude and latitude (decimal degrees)
dstart	First day to run, date in format "d/m/Y" e.g. "01/01/2016"
dfinish	Last day to run, date in format "d/m/Y" e.g. "31/12/2016"
dem	A digital elevation model used by microclima for micro-topographic effects, produced by microclima function 'get_dem' via R package 'elevatr' (internally generated via same function based on 'loc' if NA)
dem2	A digital elevation model used by microclima for meso-climate calculations, produced by microclima function 'get_dem' via R package 'elevatr' (internally generated via same function based on 'loc' if NA)
dem.res	Requested resolution of the DEM from elevatr, m
zmin	minimum elevation of DEM for terrain calculations, m (may need to be made negative if below sea level)
pixels	Number of pixels along one edge of square requested of DEM requested from elevatr, #
REFL	Soil solar reflectance, decimal %
slope	Slope in degrees (if NA, then derived from DEM with package microclima)
aspect	Aspect in degrees (0 = north) (if NA, then derived from DEM with microclima)
DEP	Soil depths at which calculations are to be made (cm), must be 10 values starting from 0, and more closely spaced near the surface
minshade	Minimum shade level to use (can be a single value or a vector of daily values) (%)
maxshade	Maximum shade level to use (can be a single value or a vector of daily values) (%)
Usrhyt	Local height (m) at which air temperature, wind speed and humidity are to be computed for organism of interest
coastal	Compute coastal effects with microclima? T (TRUE) or F (FALSE) (can take a while and may have high memory requirements depending on DEM size)
hourlydata	user input of the hourlydata matrix
dailyprecip	user input of daily rainfall
weather.elev	optional value indicating the elevation of values in 'hourlydata'. Either a numeric value, corresponding to the elevation in (m) of the location from which 'hourlydata' were obtained, or 'era5' (default, derived from Copernicus ERA5 climate reanalysis).
cad.effects	optional logical indicating whether to calculate cold air drainage effects (TRUE = Yes, slower. FALSE = No, quicker)
...	Additional arguments, see Details

Details

Parameters controlling how the model runs:

runshade = 1, Run the microclimate model twice, once for each shade level (1) or just once for the minimum shade (0)?

clearsky = 0, Run for clear skies (1) or with observed cloud cover (0)

run.gads = 1, Use the Global Aerosol Database? 1=yes (Fortran version), 2=yes (R version), 0=no

IR = 0, Clear-sky longwave radiation computed using Campbell and Norman (1998) eq. 10.10 (includes humidity) (0) or Swinbank formula (1) or from ERA5 data (2)

solonly = 0, Only run SOLRAD to get solar radiation? 1=yes, 0=no

lamb = 0, Return wavelength-specific solar radiation output?

IUV = 0, Use gamma function for scattered solar radiation? (computationally intensive)

Soil_Init = NA, initial soil temperature at each soil node, °C (if NA, will use the mean air temperature to initialise)

write_input = 0, Write csv files of final input to folder 'csv input' in working directory? 1=yes, 0=no

writecsv = 0, Make Fortran code write output as csv files? 1=yes, 0=no

windfac = 1, factor to multiply wind speed by e.g. to simulate forest

warm = 0, warming offset vector, °C (negative values mean cooling). Can supply a single value or a vector the length of the number of days to be simulated.

scenario = 0, TerraClimate climate change scenario, either 0, 2 or 4 °C warmer

terra_source = "http://thredds.northwestknowledge.net:8080/thredds/dodsC/TERRACLIMATE_ALL/data", specify location of terraclimate data, goes to the web by default

soilgrids = 0, query soilgrids.org database for soil hydraulic properties?

message = 0, allow the Fortran integrator to output warnings? (1) or not (0)

fail = nyears x 24 x 365, how many restarts of the integrator before the Fortran program quits (avoids endless loops when solutions can't be found)

spatial = 'c:/era5_data/era5', specify folder and file prefix with local ERA5 data extracted via the mcera5 package (no trailing forward slash)

save = 0, don't save forcing data (0), save the forcing data (1) or read previously saved data (2)

General additional parameters:

ERR = 1.5, Integrator error tolerance for soil temperature calculations

Refhyt = 2, Reference height (m), reference height at which air temperature, wind speed and relative humidity input data are measured

RUF = 0.004, Roughness height (m), e.g. smooth desert is 0.0003, closely mowed grass may be 0.001, bare tilled soil 0.002-0.006, current allowed range: 0.00001 (snow) - 0.02 m.

ZH = 0, heat transfer roughness height (m) for Campbell and Norman air temperature/wind speed profile (invoked if greater than 0, 0.02 * canopy height in m if unknown)

D0 = 0, zero plane displacement correction factor (m) for Campbell and Norman air temperature/wind speed profile (0.6 * canopy height in m if unknown)

Z01 = 0, Top (1st) segment roughness height(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used)

Z02 = 0, 2nd segment roughness height(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used).

ZH1 = 0, Top of (1st) segment, height above surface(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used).

ZH2 = 0, 2nd segment, height above surface(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used).

EC = 0.0167238, Eccenricity of the earth's orbit (current value 0.0167238, ranges between 0.0034 to 0.058)

SLE = 0.95, Substrate longwave IR emissivity (decimal %), typically close to 1

Thcond = 2.5, Soil minerals thermal conductivity, single value or vector of 10 specific to each depth (W/mK)

Density = 2.56, Soil minerals density, single value or vector of 10 specific to each depth (Mg/m3)

SpecHeat = 870, Soil minerals specific heat, single value or vector of 10 specific to each depth (J/kg-K)

BulkDensity = 1.3, Soil bulk density (Mg/m3), single value or vector of 10 specific to each depth

PCTWET = 0, % of ground surface area acting as a free water surface (overridden if soil moisture model is running)

rainwet = 1.5, mm of rainfall causing the ground to be 90% wet for the day

cap = 1, organic cap present on soil surface? (cap has lower conductivity - 0.2 W/mC - and higher specific heat 1920 J/kg-K)

CMH2O = 1, Precipitable cm H2O in air column, 0.1 = very dry; 1.0 = moist air conditions; 2.0 = humid, tropical conditions (note this is for the whole atmospheric profile, not just near the ground)

hori = rep(NA,24), Horizon angles (degrees), from 0 degrees azimuth (north) clockwise in 15 degree intervals

Soil moisture mode parameters:

runmoist = 1, Run soil moisture model? 1=yes, 0=no 1=yes, 0=no (note that this may cause slower runs)

PE = rep(1.1,19), Air entry potential (J/kg) (19 values descending through soil for specified soil nodes in parameter DEP and points half way between)

KS = rep(0.0037,19), Saturated conductivity, (kg s/m3) (19 values descending through soil for specified soil nodes in parameter DEP and points half way between)

BB = rep(4.5,19), Campbell's soil 'b' parameter (-) (19 values descending through soil for specified soil nodes in parameter DEP and points half way between)

BD = rep(1.3,19), Soil bulk density (Mg/m3) (19 values descending through soil for specified soil nodes in parameter DEP and points half way between)

DD = rep(2.56,19), Soil density (Mg/m3) (19 values descending through soil for specified soil nodes in parameter DEP and points half way between)

DEP and points half way between)

maxpool = 10000, Max depth for water pooling on the surface (mm), to account for runoff

rainhourly = 0, Is hourly rain input being supplied (1 = yes, 0 = no)?

rainhour = 0, Vector of hourly rainfall values - overrides daily ERA5 rain if rainhourly = 1

rainmult = 1, Rain multiplier for surface soil moisture (-), used to induce runon

rainoff = 0, Rain offset (mm), used to induce changes in rainfall from ERA5 values. Can be a single value or a vector matching the number of days to simulate. If negative values are used, rainfall will be prevented from becomming negative.

evenrain = 0, Spread daily rainfall evenly across 24hrs (1) or one event at midnight (0)

SoilMoist_Init = c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3), initial soil water content at each soil node, m3/m3

L = c(0,0,8.2,8.0,7.8,7.4,7.1,6.4,5.8,4.8,4.0,1.8,0.9,0.6,0.8,0.4,0.4,0,0)*10000, root density (m/m3), (19 values descending through soil for specified soil nodes in parameter

R1 = 0.001, root radius, m

RW = 2.5e+10, resistance per unit length of root, m3 kg-1 s-1

RL = 2e+6, resistance per unit length of leaf, m3 kg-1 s-1

PC = -1500, critical leaf water potential for stomatal closure, J kg-1

SP = 10, stability parameter for stomatal closure equation, -

IM = 1e-06, maximum allowable mass balance error, kg

MAXCOUNT = 500, maximum iterations for mass balance, -

LAI = 0.1, leaf area index (can be a single value or a vector of daily values), used to partition traspiration/evaporation from PET in soil moisture model

microclima.LAI = 0, leaf area index, used by package microclima for radiation calcs

LOR = 1, leaf orientation for package microclima radiation calcs

Snow mode parameters:

snowmodel = 1, run the snow model 1=yes, 0=no (note that this may cause slower runs)

snowtemp = 1.5, Temperature (°C) at which precipitation falls as snow

snowdens = 0.375, snow density (mg/m3), overridden by densfun

densfun = c(0.5979, 0.2178, 0.001, 0.0038), slope and intercept of model of snow density as a linear function of snowpack age if first two values are nonzero, and following the exponential function of Sturm et al. 2010 J. of Hydromet. 11:1380-1394 if all values are non-zero; if it is c(0,0,0,0) then fixed density used

snowmelt = 1, proportion of calculated snowmelt that doesn't refreeze

undercatch = 1, undercatch multipier for converting rainfall to snow

rainmelt = 0.0125, paramter in equation that melts snow with rainfall as a function of air temp

snowcond = 0, effective snow thermal conductivity W/mC (if zero, uses inbuilt function of density)

intercept = max(maxshade) / 100 * 0.3, snow interception fraction for when there's shade (0-1)

grasshade = 0, if 1, means shade is removed when snow is present, because shade is cast by grass/low shrubs

Intertidal mode parameters:

shore Include tide effects? If 1, the matrix tides is used to specify tide presence, sea water temperature and presence of wavesplash

tides = matrix(data = 0, nrow = length(seq(as.POSIXct(dstart, format = '

Outputs:

ndays - number of days for which predictions are made

longlat - longitude and latitude for which simulation was run (decimal degrees)

dates - vector of dates (POSIXct, UTC)

nyears - number of years for which predictions are made

RAINFALL - vector of daily rainfall (mm)

elev - elevation at point of simulation (m)

minshade - minimum shade for simulation (%)

maxshade - maximum shade for simulation (%)

dem - digital elevation model obtained via 'get_dev' using package 'elevatr' (m)

DEP - vector of depths used (cm)

SLOPE - slope at point of simulation (%)

ASPECT - aspect at point of simulation (°, 0 is north)

HORIZON - horizon angles at point of simulation (°)

metout/shadmet variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3 TALOC - air temperature (°C) at local height (specified by 'Usrhyt' variable)

• 4 TAREF - air temperature (°C) at reference height (specified by 'Refhyt', 2m default)

• 5 RHLOC - relative humidity (%) at local height (specified by 'Usrhyt' variable)

• 6 RH - relative humidity (%) at reference height (specified by 'Refhyt', 2m default)

• 7 VLOC - wind speed (m/s) at local height (specified by 'Usrhyt' variable)

• 8 VREF - wind speed (m/s) at reference height (specified by 'Refhyt', 2m default)

• 9 SNOWMELT - snowmelt (mm)

• 10 POOLDEP - water pooling on surface (mm)

• 11 PCTWET - soil surface wetness (%)

• 12 ZEN - zenith angle of sun (degrees - 90 = below the horizon)

• 13 SOLR - solar radiation (W/m2) (unshaded, adjusted for slope, aspect and horizon angle)

• 14 TSKYC - sky radiant temperature (°C)

• 15 DEW - dew fall (mm / h)

• 16 FROST - frost (mm / h)

• 17 SNOWFALL - snow predicted to have fallen (cm)

• 18 SNOWDEP - predicted snow depth (cm)

• 19 SNOWDENS - snow density (g/cm3)

soil and shadsoil variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3-12 D0cm ... - soil temperature (°C) at each of the 10 specified depths

if soil moisture model is run i.e. parameter runmoist = 1

soilmoist and shadmoist variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3-12 WC0cm ... - soil moisture (m3/m3) at each of the 10 specified depths

soilpot and shadpot variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3-12 PT0cm ... - soil water potential (J/kg = kPa = bar/100) at each of the 10 specified depths

humid and shadhumid variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3-12 RH0cm ... - soil relative humidity (decimal %), at each of the 10 specified depths

plant and shadplant variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3 TRANS - plant transpiration rate (g/m2/h)

• 4 LEAFPOT - leaf water potential (J/kg = kPa = bar/100)

• 5-14 RPOT0cm ... - root water potential (J/kg = kPa = bar/100), at each of the 10 specified depths

if snow model is run i.e. parameter snowmodel = 1
sunsnow and shdsnow variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3-10 SN1 ... - snow temperature (°C), at each of the potential 8 snow layers (layer 8 is always the bottom - need metout$SNOWDEP to interpret which depth in the snow a given layer represents)

if wavelength-specific solar output is selected i.e. parameter lamb = 1
solar output variables drlam (direct solar), drrlam (direct Rayleigh solar) and srlam (scattered solar) variables:

• 1 DOY - day-of-year

• 2 TIME - time of day (mins)

• 3-113 290, ..., 4000 - irradiance (W/(m2 nm)) at each of 111 wavelengths from 290 to 4000 nm

Value

metout The above ground micrometeorological conditions under the minimum specified shade

shadmet The above ground micrometeorological conditions under the maximum specified shade

soil Hourly predictions of the soil temperatures under the minimum specified shade

shadsoil Hourly predictions of the soil temperatures under the maximum specified shade

soilmoist Hourly predictions of the soil moisture under the minimum specified shade

shadmoist Hourly predictions of the soil moisture under the maximum specified shade

soilpot Hourly predictions of the soil water potential under the minimum specified shade

shadpot Hourly predictions of the soil water potential under the maximum specified shade

humid Hourly predictions of the soil humidity under the minimum specified shade

shadhumid Hourly predictions of the soil humidity under the maximum specified shade

plant Hourly predictions of plant transpiration, leaf water potential and root water potential under the minimum specified shade

shadplant Hourly predictions of plant transpiration, leaf water potential and root water potential under the maximum specified shade

sunsnow Hourly predictions of snow temperature under the minimum specified shade

shadsnow Hourly predictions snow temperature under the maximum specified shade

Examples

library(NicheMapR)
library(ecmwfr)
library(mcera5)
library(lubridate)
library(dplyr)
library(tidync)

# get ERA5 data with package mcera5 (just do once for region and time of interest)

# assign your credentials (register here: https://cds.climate.copernicus.eu/user/register)
uid <- "$$$$$$"
cds_api_key <- "$$$$$$$$-$$$$-$$$$-$$$$-$$$$$$$$$$$$"

ecmwfr::wf_set_key(user = uid, key = cds_api_key, service = "cds")

# bounding coordinates (in WGS84 / EPSG:4326)
xmn <- 130
xmx <- 132
ymn <- -26
ymx <- -24

# temporal extent
st_time <- lubridate::ymd("2010:07:01")
en_time <- lubridate::ymd("2011:12:31")

# filename and location for downloaded .nc files
file_prefix <- "era5"
op <- "C:/Spatial_Data/"

# build a request (covering multiple years)
req <- build_era5_request(xmin = xmn, xmax = xmx,
                          ymin = ymn, ymax = ymx,
                          start_time = st_time,
                          end_time = en_time,
                          outfile_name = file_prefix)
str(req)
request_era5(request = req, uid = uid, out_path = op)

# run micro_era5 for a location (make sure it's within the bounds of your .nc files)

dstart <- "01/01/2011"
dfinish <- "31/12/2011"
loc <- c(131, -25) # somewhere in the middle of Australia
micro<-micro_era5(loc = loc, dstart = dstart, dfinish = dfinish, spatial = 'c:/Spatial_Data/era5')

metout<-as.data.frame(micro$metout) # above ground microclimatic conditions, min shade
soil<-as.data.frame(micro$soil) # soil temperatures, minimum shade
soilmoist<-as.data.frame(micro$soilmoist) # soil temperatures, minimum shade

# append dates
tzone<-paste("Etc/GMT+",0,sep="")
dates<-seq(as.POSIXct(dstart, format="%d/%m/%Y",tz=tzone)-3600*12, as.POSIXct(dfinish, format="%d/%m/%Y",tz=tzone)+3600*11, by="hours")

metout <- cbind(dates,metout)
soil <- cbind(dates,soil)
soilmoist <- cbind(dates, soilmoist)

# plotting above-ground conditions in minimum shade
with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)"
, type = "l",main=paste("air and sky temperature",sep=""), ylim = c(-20, 60))})
with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)"
, type = "l",lty=2,col='blue')})
with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)"
,  type = "l",col='light blue',main=paste("sky temperature",sep=""))})
with(metout,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (%)"
, type = "l",ylim=c(0,100),main=paste("humidity",sep=""))})
with(metout,{points(RH ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (%)"
, type = "l",col='blue',lty=2,ylim=c(0,100))})
with(metout,{plot(VREF ~ dates,xlab = "Date and Time", ylab = "Wind Speed (m/s)"
,  type = "l",main="wind speed",ylim = c(0, 15))})
with(metout,{points(VLOC ~ dates,xlab = "Date and Time", ylab = "Wind Speed (m/s)"
,  type = "l",lty=2,col='blue')})
with(metout,{plot(SOLR ~ dates,xlab = "Date and Time", ylab = "Solar Radiation (W/m2)"
,  type = "l",main="solar radiation")})
with(metout,{plot(SNOWDEP ~ dates,xlab = "Date and Time", ylab = "Snow Depth (cm)"
,  type = "l",main="snow depth")})

# plotting soil temperature
for(i in 1:10){
 if(i==1){
   plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)"
   ,col=i,type = "l",main=paste("soil temperature",sep=""))
 }else{
   points(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature
    (°C)",col=i,type = "l")
 }
}

# plotting soil moisture
for(i in 1:10){
 if(i==1){
   plot(soilmoist[,i+3]*100~soilmoist[,1],xlab = "Date and Time", ylab = "Soil Moisture (% volumetric)"
   ,col=i,type = "l",main=paste("soil moisture",sep=""))
 }else{
   points(soilmoist[,i+3]*100~soilmoist[,1],xlab = "Date and Time", ylab = "Soil Moisture
    (%)",col=i,type = "l")
 }
}

mrke/NicheMapR documentation built on May 5, 2024, 1:13 a.m.