Microclima: An r package for modelling meso‐ and microclimate

Climate is of fundamental importance to the ecology and evolution of all organisms. However, studies of climate–organism interactions usually rely on climate variables interpolated from widely spaced measurements or modelled at coarse resolution, whereas the conditions experienced by many organisms vary over scales from millimetres to metres. To help bridge this mismatch in scale, we present models of the mechanistic processes that govern fine‐scale variation in near‐ground air temperature. The models are flexible (enabling application to a wide variety of locations and contexts), can be run using freely available data and are provided as an R package. We apply a mesoclimate model to the Lizard Peninsula in Cornwall to provide hourly estimates of air temperature at resolution of 100 m for the period Jan‐Dec 2010. A microclimate model is then applied to a 1 km2 region of the Lizard Peninsula, Caerthillean Valley (49.969°N, 5.215°W), to provide hourly estimates of near‐ground air temperature at resolution of 1 m2 during May 2010. Our models reveal substantial spatial variation in near‐ground temperatures, driven principally by variation in topography and, at the microscale, by vegetation structure. At the meso‐scale, hours of exposure to air temperatures at 1 m height in excess of 25°C ranged from 23 to 158 hr, despite this temperature never being recorded by the weather station within the study area during the study period. At the micro‐scale, steep south‐facing slopes with minimal vegetation cover experienced temperatures in excess of 40°C. The microclima package is flexible and efficient and provides an accurate means of modelling fine‐scale variation in temperature. We also provide functions that facilitate users to obtain and process a variety of freely available datasets needed to drive the model.

of kilometres. In contrast, the conditions experienced by many organisms vary over scales from millimetres to metres (Potter, Woods, & Pincebourde, 2013). This spatial mismatch is bridged implicitly in many models by assuming that grid-cell average climatic variables are statistically meaningful predictors of ecological responses (Bennie, Wilson, Maclean, & Suggitt, 2014). Statistical associations between organism and coarse-gridded climate data are therefore widely used, and have shown themselves be powerful predictive tools (Guisan & Thuiller, 2005). However, to investigate mechanistic links between climate and physiology, the effects of short-term variability and the role of microclimates in buffering ecological change, fine-resolution data are required. Thus, much ecological and evolutionary research is still hampered by an inability to model climate at fine resolution (Potter et al., 2013;Suggitt et al., 2017).
Despite the tendency for ecologists to use coarse-resolution climate data, studies of microclimate have a long history and many of the processes were well understood more than 30 years ago (Campbell, 1986;Geiger, 1927;Hay, 1979). However, many of these early studies drew on field measurements and studied aspects of microclimate at single locations. Ecologists often require data over larger spatial extents, and gridded climate data are particular useful (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). Recent advances in remote-sensing and the growing availability of very fine-resolution remotely derived datasets create a timely opportunity to present methods and models capable of generating gridded climate datasets at fine-resolution.
To date, several approaches to downscaling from coarse-gridded to fine-scale microclimate data have been used by ecologists.
Dynamical downscaling, through the use of regional climate models that apply the full physics of global climate models at a fine-scale (Murphy, 2000), have the advantage that they can generate internally consistent data for variables and represent synoptic systems.
However, due to high computing requirements they are rarely a practical solution for producing data at resolutions below five km.
Physically based boundary-layer models of atmospheric processes at finer scales (e.g. down to 1 m resolution) are usually limited in application to a small extent and highly simplified landscapes. Land surface models (e.g. JULES, the UK land surface simulator) apply physical equations to solve the energy and water balance at a point, or across a grid, and in doing so predict key ecological variables such as near-surface temperature and humidity. However, while land surface models incorporate vertical processes such as radiative heating of the surface and canopy shading, they do not incorporate mesoscale processes such as variation in wind speed, elevational lapse rates or lake/ocean effects. While land surface models have been adapted for use in an ecological context (Bennie et al., 2010), most physically based models are primarily designed for meteorological or hydrological applications. A notable exception is the NicheMapR package in r (Kearney & Porter, 2017), which is explicitly designed to mechanistically model the energy and mass budgets of organisms and their microclimate (including soil and snow), and has been widely tested (see e.g. Kearney et al., 2014). Finally, GIS-based statistical downscaling techniques apply empirical corrections (usually based on slope, aspect and elevation) to map climatic variables, and have been used in several studies to produce fine-resolution maps for species distribution modelling (e.g. Milling et al., 2018).
The models and r package described here are not intended to replace physically based regional climate models, land surface schemes or mechanistic approaches to the energy budget of organisms and their environment such as NicheMapR. Of note, however, many of required to model near-ground temperature are similar to those required for modelling the heat budget organisms, and the functions in microclima may be of use in so doing. However, our primary intention is to bridge the gap between the landscape and local-scale processes that cause spatial variation in temperature and can be modelled using fine-resolution Digital Terrain Models (DTMs) and point-based models to determine the energy balance (Table 1). We develop a flexible hybrid physically-and empirically based approach in which the spatial patterns of physical factors directly influencing the near-ground temperatures at a point are calculated, and the relative influence of these factors within a given landscape or region can be fitted to data by empirically derived parameters. This hybrid approach to mapping climatic variables at a fine scale is suitable for many ecological applications, avoiding the complexity and computational costs of attempting to fully resolve the physics of atmospheric processes at high resolution, but retaining much of the generality of a physically based model. The models are designed to be flexible, enabling application in a wide variety of circumstances, though their modular design is such that easy development of improvements for application in specific circumstances is also possible. The models can also be easily applied TA B L E 1 Summary of modelling approaches used for microclimate research

| Radiation
The net radiation flux is determined by the balance of incoming shortwave radiation and emitted longwave radiation, with the former portioned between direct (R dir ) and diffuse (R dif ) components.

Shortwave radiation is modified by topography and vegetation cover
and downscaled using function shortwaveveg. Topography determines whether a given location is shaded and also the angle at which the sunlight strikes the surface. Vegetation attenuates radiation as it passes though the canopy.
From Bennie et al. (2008), the direct radiation flux on an inclined surface is given by: where R beam is the direct beam radiation flux on a surface perpendicular to the beam, Z is the solar zenith, S is the angle of the slope of the surface, Ω S is the solar azimuth, Ω is the slope aspect, (1) A is the solar altitude and H is the local horizon angle in the direction of the sun. Z, A and Ω S can be readily determined for a given time and geographical position and the slope and aspect of a surface and local horizon angles from digital elevation data.
From Hay and McKay (1985), the diffuse radiation flux can be partitioned into that which is isotropically distributed (R * dif ), that which exhibits anistropic properties (R ′ dif ) and that which is reflected where α S is the mean albedo of the surrounding surface and S* is the mean slope of the adjacent surfaces. The relative partitioning of radiation depends on an "anisotropy index" (k) given by: where R 0 is the extraterrestrial radiation flux (~4.87 MJ m −2 hr −1 ) and s is a correction for the proportion of sky, calculated using function skyviewtopo, as follows: where H is the mean horizon angle.
The transmission of radiation by vegetation is described using an equation similar to Beer's Law (Campbell, 1986): where R veg is the flux density of radiation absorbed by the ground below leaf-area index (L AI ), α g is the albedo of the ground below the canopy, K′ and K* are the isotropic and anisotropic coefficients of the canopy and s veg is an adjustment applied if the sky view above the canopy is partially obscured (see later). K′ is a function of solar inclination and leaf distribution character of the canopy. From Campbell (1986), a broad range of leaf types can be represented by an ellipsoidal distribution, and the extinction coefficient can thus be expressed as follows: Here, x is determined by canopy architecture and is the ratio of vertical to horizontal projections of a representative volume of foliage, and in our model is estimated allometrically from vegetation height using function leaf_geometry (Appendix S2). The extinction coefficient for isotropic component of radiation (K*) can be obtained by integrating over the portion of the hemisphere in view. For computational efficiency, the integral can be closely approximated by Equation 3, with A (in degrees) substituted by a parameter A* which, for a given values of x, can be derived from L AI as follows: where p 1 and p 2 are coefficients unique to each x (Supporting Information Table S3). If the sky view above the canopy is partially obscured, then the integral is between the limits determined by H and the sky view correction factor (s veg ) is applied. In function skyviewveg, this integral is approximated by Equation 2, with H replaced by H*: Here, p 3 and p 4 are parameters unique to each value of x (Supporting Information Table S4).
Following Konzelmann et al. (1994), and assuming that differences in R lw caused by difference between T and T 0 are small, the net flux of longwave radiation under vegetated canopies (R lw ), calculated using function longwaveveg, can be approximated as follows: where σ is the Stefan and Boltzmann constant (2.043 × 10 −10 MJ m −2 hr −1 ), R lwg is radiation emitted back from the atmospheric that passes through gaps the canopy, R lwe is radiation scattered downwards from leaves, R lwc is radiation emitted by the canopy and T is temperature in Kelvin.
The flux of radiation that passes through gaps in the canopy is given by: R lsky is radiation scattered back from the atmosphere which, assuming that differences in R lw caused by differences between T and T 0 are small, can be calculated as follows: Here, ε is the emissivity of the atmosphere, which can be determined as follows (Klok & Oerlemans, 2002): where n is fractional cloud cover and e a is vapour pressure in kPa.
From Zhao and Qualls (2006), the flux of radiation scattered downward through leaf reflection is given by: where α c is the albedo of the canopy and r is the fraction of downward radiation scattered upwards, estimated as: R lwc is given by:

| Wind speed
Wind speeds are affected by local terrain, and to account for this, function wind-height implements the logarithmic wind-height profile assumed by Allen, Pereira, Raes, and Smith (1998), and function windcoef applies the shelter coefficient described by Ryan (1977), as follows: where u 1 is local wind speed at 1 m above the ground, u 10 is the wind speed at 10 m height and H w is the tangent of the horizon angle upwind at 1 m above the ground.

| Mesoclimate model
The mesoclimate model provides estimates of air temperature and ignores the effects of radiation transmissions though canopies and variation in ground surface albedo, as these are accounted for in the microclimate model. Differences between local temperatures (T 1 ) and reference air temperature (T) are derived as a function of coastal, cold-air drainage and elevation effects and also the effects of meso-scale topography on the radiation flux, as in Equation 1: here, ΔT E is the difference in temperature due to elevation, ΔT C is the difference in temperature due to coastal effects and ΔT K is the difference in temperature due to cold-air drainage.

| Elevation
Differences in temperature due to elevation are calculated as follows: where Δz is the difference in elevation (m) between the locations of T and T 1 and Γ m is the lapse rate, calculated using function lapserate, as follows (Hess, 1959): where g is gravitational acceleration (9.8076 m/s), L v is the latent heat of vaporisation (2,501,000 J/kg), Q is the gas constant for dry air (287 J kg −1 K −1 ), c pd is the specific heat of dry air at constant pressure (1,003.5 J kg −1 K −1 ), T is the reference temperature and r v is the mixing ratio of water vapour given by: where P is atmospheric pressure (Pa).

| Coastal effects
Coastal effects are derived using function coastalTps, which uses thin-plate spline interpolation with three covariates to derive finer resolution temperature estimates for each time step from coarsegridded reference temperature data. The three covariates are: differences between sea and reference land temperature, coastal exposure in an upwind direction and coastal exposure irrespective of direction. Upwind exposure is calculated as the inverse-square distance weighted proportion of sea upwind of each location and general exposure by numerically integrating this ratio at fixed intervals over the full 360°.

| Cold-air drainage
ΔT K is modelled as follows: where I C is a binary variable, conditional on time of day, wind speed and emissivity, as cold-air drainage typically occurs at night or shortly after, and in calm, still conditions (Barr & Orgill, 1989).
The function cadconditions, used for calculating I C allows the user to specify these conditions. Δz m is the elevation difference in metres of a given location and the highest point of a drainage basin, and F is accumulated flow expressed as a proportion of the maximum in each basin, and calculated using function flowacc.
Quantification of F and Δz m requires drainage basins to be delineated, using function basindelin.

| Data
To calibrate and run the models, the following high spatial, low temporal resolution datasets are needed (summarised in (5) u 1 = 0.635u 10 1 − arctan 0.17H w 1.65

| Vegetation characteristics
Following for example Carlson and Ripley (1997), we estimated leaf-area index from the normalised difference vegetation index (NDVI), using visual and colour-infrared aerial imagery obtained

| Albedo
We derived albedo from the same visual and colour-infrared aerial imagery, adjusting values for brightness and contrast using MODIS data obtained from USGS Land Processes Distributed Archive Centre (Appendix S2 and functions albedo and albedo _ adjust).

| Cloud cover and shortwave radiation
We used 0.05º gridded satellite-derived estimates of cloud cover, and direct and diffuse radiation (Posselt et al., 2012). Radiation data are available hourly, though missing values and those within an hour either side of sunrise and sunset, for which satellite estimates are unreliable (Posselt et al., 2012), were imputed (Appendix S2). Cloud cover is available at ~15 min intervals and each grid cell is assigned a value of "full," "partial" or "unobscured." Fractional cloud cover was calculated by calculating the mean in each hour, making the assumption that partial cloud cover equates to fractional value of 0.5.

| Surface pressure and wind data
We obtained six-hourly surface pressure and wind data from the National Center for Environmental Prediction (Kalnay et al., 1996).
These data are available at a grid cell resolution of 2.5°, and the values for the grid cell corresponding to our study area were extracted.
Values were then interpolated to hourly data using a cubic-spline.

| Humidity and temperature
Daily specific humidity data, mean daily near-surface air temperature, and daily temperature ranges, available at a 1 km grid resolution, were obtained from the Centre for Ecology and Hydrology Climate (Robinson et al., 2015). Hourly specific humidity data were derived by interpolation using a cubic-spline. To derive hourly temperature data, we implemented a more complex interpolation algorithm, whereby diurnal patterns and variation in cloud cover and radiation are accounted for (Appendix S2 and function hourlytemp).

| Sea-surface temperature
We obtained 1 degree gridded datasets of monthly sea ice and sea-surface temperatures, available as a global dataset from the Met Office Hadley Centre (Rayner et al., 1996), and extracted data for the grid cell corresponding to our study area. We obtained hourly values using cubic-spline interpolation, assuming that the mean value for each month corresponded to the midpoint of that month. Due to the high volume and specific heat capacity of water, sea-surface temperatures undergo only minor high frequency fluctuations, so simple interpolation was deemed reasonable.

| Model fitting
Prior to fitting the mesoclimate model we accounted for cold-air drainage, elevation and coastal effects. To calculate elevation effects, we first removed the fixed lapse rate applied to the temperature data and then applied our variable one.
To fit our mesoclimate model, 56 iButton thermachrons were deployed across the Lizard Peninsula between 1 March and 31 December 2010, and set to record temperatures at hourly intervals. Loggers were placed to capture the full spatial gradients in the main determinants of climate to minimise extrapolation errors, and provided 137,218 measurements of temperature.
Loggers were attached to a wooden pole 1 m above the ground.
To fit the microclimate model, 55 iButton thermachrons were deployed in Caerthillean Valley on the Lizard Peninsula (49.9687ºS, 5.2142ºW), from 10 May to 31 May 2010. Loggers were set to record temperatures at 10 min intervals, and the mean temperature in each hour used to calibrate the model. 27,530 hourly temperature measurements were obtained. The valley is a coastal grassland with complex topography, enabling temperatures to be recorded across a wide range of slopes and aspects and in vegetation of varying height. Loggers were attached to a short wooden stake 5 cm above the ground. In both instances, loggers were orientated north, and shielded from direct sunlight using a white plastic screen. Data from half the loggers were used for calibration and from the rest for testing.
Temperature anomalies (T 0 − T) were modelled using standard linear regression as a function of the following sets of terms: Here, u f is a factor that allows the relationship with net radiation to vary with wind speed, set at 0 when wind speeds at 1 m are (7) T 0 − T = β 1 +β 1 R net +β 1 u f +β 1 R net u f +β 1 ΔT K +ε i below 3.66 m/s, and one when above (β 4 is assumed to be negative), with this threshold established by iteratively trying out different thresholds, and selecting that which yielded the best fit. The terms β 1,…,5 are coefficients estimated by linear regression and ε the error associated with each term i. Other terms have already been defined.
The terms are listed in anticipated descending order of importance.
We first assessed whether including each set of terms improved model performance by computing the Akaike information criterion (AIC) and then estimated coefficients associated with each term using standard linear regression. To reduce the effects of temporal autocorrelation, we randomly selected 2,000 of the temperature measurements and repeated the analyses 9,999 times, computing AICs and coefficient estimates for each model run. Median model coefficient estimates were then used to drive our model. The microclimate model was fitted in the same way, except that here (T 0 − T) is the difference in near ground temperatures at the output of the mesoclimate model and the value of u f that yielded the best fit was 0.398 m/s. The function fitmicro implements the method described above, though also includes the option to use all data for fitting. taken to run the model on a 1,000 × 1,000 pixel dataset took 0.25 s for one time-step, equating to just 36 min for a year (though additional time is required to write model outputs to disk).

| RE SULTS
In all 9,999 model simulations, both sets of models performed best when all terms were included. This mesoclimate model explained 90.8% of the variation in local temperature anomalies and 96.2% of the variation in total temperature, with a mean absolute error (MAE) of 0.97°C and root mean square error (RMSE) of 1.23°C (Figure 2a,c). The microclimate model explained 78.7% of the variation in local temperature anomalies and 90.9% of the variation in total temperature, with a MAE of 1.25°C and RMSE of 1.61°C (Figure 2a,c). Model coefficients for both the meso-and microclimate model are shown in Table 2.
At the meso-scale, there was relatively little spatial variation in mean temperature, which in 2010 ranged from 8.6°C to 10.0°C ( Figure 3a). The warmest temperatures were in sheltered low-lying coastal valleys, particularly on south-facing slopes. Minimum temperatures ranged from −6.7°C to −5.3°C, being coldest at higher elevations and inland (Figure 3b). Maximum temperatures ranged from 27.2°C to 31.2°C, and were highest on low-lying south-facing slopes (Figure 3c).
There were larger differences in bioclimate variables. Accumulated degree hours varied from 8,446 to 16,008 (Figure 3d), hours of exposure F I G U R E 2 Observed and predicted temperatures. In (a) temperatures recorded by iButtons place 1 m above the ground are compared to outputs obtained from the mesoclimate model, and in (b) temperatures recorded by iButtons 5 cm above-ground level are compared to the outputs of the microclimate model. In  (Figure 4c). There were also large differences in bioclimate variables. Growing-degree hours varied from 1,644 to 4,223 (Figure 4d), mean diurnal temperature variation from 11.1°C to 21.3°C (Figure 4e) and hours of exposure to temperatures in excess of 30°C from 0 to 53 ( Figure S4).

| D ISCUSS I ON
The main aim of this study was to present general methods for modelling micro-and mesoclimate that can be readily applied to determine the range in near-ground air temperatures experienced by organisms across any landscape or region. While the models accurately predict temperatures at locations other than those used for model calibration, their transferability to different sites altogether has yet to be tested and, although calibration and testing were performed under a wide range of climatic conditions, there may be errors associated with extrapolating the model beyond the conditions used for calibration. However, an important characteristic of our models is that the spatial patterns of variables are based on the underlying physics of heat budgets and airflow rather than on spatial interpolation, and while recalibration or the incorporation of other macro-to microscale processes may be necessary at some locations, the physical laws governing these processes are universal.
Overall, the predictive power of our models compare well with other more location-specific models (Aalto, Riihimäki, Meineri, Hylander, & Luoto, 2017;Pike, Pepin, & Schaefer, 2013), and build TA B L E 2 Median, mean (±1 SD) model coefficients associated with meso-and microclimate model  Bennie et al., 2008) or by incorporating mesoclimatic processes (cf. Kearney & Porter, 2017). Nonetheless, some aspects of the model remain poorly developed, in part due to the limited extent over which it has been tested, and hence, the range of conditions that influence climatic processes within our study area. Key among these is the effects of latent heat flux on temperatures, which can be particularly important in cold environments, where snow freeze-thaw is frequent (Weller & Holmgren, 1974), or under drought conditions when soil temperatures may heat up by more than predicted (Hunt, Kelliher, McSeveny, & Byers, 2002). In contrast to other models (e.g. Kearney & Porter, 2017), heat exchange between the soil and near-ground air layer and heat storage in the soil are also unaccounted for, and may result in delayed effects of radiation on near-ground temperatures. Environmental lapse rates are also rather crudely handled by our model; for our study area, this does not cause large errors due to the limited elevation range, but further development and testing may be necessary for applications in mountainous regions. A further limitation is that our model does not presently account for seasonal variation in albedo, which in temperate regions can be significant due to leaf-loss in winter, and in Arctic regions may be influenced strongly by snow cover (Aalto et al., 2017;Weller & Holmgren, 1974). Vegetation structure is also rather simplistically determined from aerial imagery. Better three-dimensional assessment of seasonal variation in vegetation structure, made possible through fullwaveform laser-scanning, for example (Wagner, Hollaus, Briese, & Ducic, 2008), represents one of the best opportunities for further development. These limitations aside, our models provide accurate physically based predictions of the effects of topography and vegetation on local-scale climate at the landscape scale.
At both micro-and meso-scales, slope and aspect are the principal determinants of spatial variation in maximum temperatures, with the warmest temperatures on steep south-facing slopes. However, at the micro-scale, where surface albedo and vegetation structure are also accounted for, these also exert a strong influence, with temperatures highest on dark surfaces with sparse vegetation cover. This is to be expected given the overriding importance of net solar radiation on temperature (Geiger, 1927 (Maclean, Suggitt, Wilson, Duffy, & Bennie, 2017). This is largely due to the particularly cold winter that affected much of north and north-western Europe, caused by record persistence of the negative phase of the North-Atlantic Oscillation (Cattiaux et al., 2010). The total number of frost-hours (<0°C) recorded at Culdrose was the greatest on record, more than eight times higher than the 1977-2016 median.
This is reflected in spatial patterns of frost exposure across the study region, which even in sheltered valleys is relatively high, despite being frost frost-free in many years (Maclean et al., 2017). Rather uncharacteristically, the maximum recorded tem- Biological responses to climate change within our study region are influenced strongly by fine-scale spatial and temporal variation (Maclean, Hopkins, Bennie, Lawson, & Wilson, 2015). Consequently, predictions of the responses of species to climate change will need to account for the spatial variation in microclimate at resolution smaller than most available climate data, and the dynamics of microclimate at a temporal resolution smaller than long-term climatic means. More generally, the study of relationships between species and climate is currently hampered by the coarse resolution at which climate is currently modelled (Bramer et al., 2018;Potter et al., 2013;Suggitt et al., 2018). This study is intended to demonstrate the importance of fine-scale variation in temperature and to show that this variation can be modelled.

ACK N OWLED G EM ENTS
We thank Michael Kearney and two anonymous referees for useful comments on the manuscript, Phillipa Gillingham, Robin Curtis and Isobel Bramer for help with fieldwork and Robert Wilson for conceptual advice. We also thank Katie Thurston for proof-reading, though any remaining mistakes are the authors' responsibility. This research was part-funded by the European Social Fund Project 09099NCO5 and NERC NE/P016790/1.

DATA ACCE SS I B I LIT Y
The data used in this study are included with the r package, available at https://doi.org/10.5281/zenodo.1411517 (Maclean, 2018).