The influence of system heterogeneity on peat-surface temperature dynamics

Temperatures at the soil–atmosphere interface influence ecosystem function by driving nonlinear terrestrial biogeochemical, ecohydrological, and micrometeorological processes. Whilst climate, soil and vegetation controls on spatially average ecosystem temperatures are recognised, how interacting and heterogeneous ecosystem layers create spatio-temporal complex thermal ecosystems has not been determined. Such thermal hot spots and hot moments may underpin the capability of ecosystems to support biological and biogeochemical diversity and control the likelihood of tipping points in system-regulating feedbacks being locally exceeded. This is of notable importance in peatlands, where soil temperatures control the storage of their associated globally important carbon stocks. Here, through the application of high spatio-temporal resolution surface temperature data and peat thermal modelling, we assess the impact of system heterogeneity (spatio-temporal impact of the following system layers: tree, shrubs, microtopography, groundcover species and sub-surface ice) on surface temperature regimes. We show (a) that peat-surface thermal hotspot intensity and longevity is linked to system heterogeneity and (b) that not all system layers have an equal influence over the peat-surface thermal regime and extreme temperatures; thermal heterogeneity increases up to a maximum of five layers of heterogeneity and decreases thereafter. The results crucially demonstrate that such changes in the spatio-temporal thermal dynamics and extremes may occur without significant changes in median temperatures. This is important to the conceptual understanding of peatland responses and ecosystem resilience to disturbance. It emphasises the need to determine the potential for transitions in magnitude, longevity and locality of small-scale thermal extremes to induce functional transitions that propagate through given ecosystems, and to characterise the impact of such small-scale spatio-temporal complexity on ecosystem scale biogeochemical and ecohydrological function.


Introduction
The soil-surface is a critical interface, controlling mass and energy exchange between terrestrial ecosystems and the atmosphere. Dynamic, spatially complex soil-surface temperatures regulate the rates of such interface exchanges, governing processes such as carbon storage and release (Kirschbaum 1995 (Frei et al 2012). Surface temperatures also influence processes at depth by driving the thermal regime of the soil profile (Wullschleger et al 1991, Kettridge et al 2013. This control of system processes by surface temperature is of notable importance within peatland ecosystems where complex peat-surfaces, which are impacted extensively by varying disturbance regimes (Turetsky et al 2002), act to regulate deeper stores of organic carbon (Belyea andClymo 2001, Belyea andBaird 2006) that account for one-third of the global soil carbon pool (Turunen et al 2002).
Despite a past focus on point measurements (Kellner 2001, Kellner and Halldin 2002, Lafleur et al 2005, Kettridge and Baird 2007, Mckenzie et al 2007a, Kettridge et al 2012, peatland surface temperatures are rarely uniform. Peat-surface temperatures can vary by over 25 • C at the decimetre scale in undisturbed, forested conditions (Leonard et al 2018). This complex thermal signature emerges from interactions between the meteorological regime, energy interception by vegetation (Hardy et al 2004, Leonard et al 2018, and the organisation of geomorphological, hydrological (Al-Kayssi 2002, Folwell et al 2015, thermal (Peters-Lidard et al 1998) and aerodynamic (Green et al 1984) properties. The total structural complexity of a peatland atmosphere interface may be conceptualised as vertical (variability in structure across a vertical space) and horizontal heterogeneity (variability across each layer found in the vertical space (Kelliher et al 1993, Filotas et al 2014). For example, the vertical layers of a peatland with tree cover comprises of: (a) tree layer, (b) a vascular sub canopy of ericaceous shrubs, (c) a moss groundcover (of Sphagnum mosses, true mosses, lichen and patches of bare peat (d) microtopography and sometimes (e) ice layers. Each of these vertical layers varies in space within each layers' horizontal plane. The peat surface temperature is the sum of complicated small-scale, spatio-temporal interactions across and between its horizontal and vertical layers. For example, vascular species such as trees and shrubs influence the energy available at lower system layers. The distinctive micro-topography of hummocks and hollows found in peatlands (Sjörs 1961, Foster et al 1983 can induce surface temperature variations of up to 13 • C between north-facing and south-facing hummock slopes less than a meter apart. This variability is primarily a result of sun angle and direct shortwave radiation input (Kettridge and Baird 2010), but micro-topography is also intrinsically linked to water table depths and soil moisture which further modify peat temperatures (Kettridge and Baird 2008). Peat ground-cover also has an important influence over peat-surface (Kettridge and Baird 2010) temperatures (Soudzilovskaia et al 2013), with 20 • C temperature variations observed between neighbouring patches of Pleaurozium schreberi and Sphagnum fuscum (Stoy et al 2012). Processes inducing such temperature differences between ground-cover types include large differences in surface albedo (Stoy et al 2012) and evaporative cooling (Brown et al 2010, Blok et al 2011. Vegetation may control evaporative cooling through functional differences in their ability to meet evaporative demand. For example, Sphagnum species can access water from depth by means of capillary rise, unlike feather mosses, leading to differences in surface resistance to evaporation between species (Mccarter andPrice 2014, Leonard et al 2017). Subsurface heterogeneity may also influence surface thermal regimes when seasonally frozen, uneven ice thawing results in the presence of frozen peat and frost-free peat patches well into the growing season (Petrone et al 2008). With the spatial variations in peatland surface temperatures similar in magnitude to diurnal fluctuations (Leonard et al 2018), understanding the variation in the thermal response of peatlands to climatic shifts or more localised disturbance is likely as important as determining the average system behaviour.
Whilst both the complex spatio-temporal distribution in peat surface temperatures and the processes that induce them have been highlighted above, the magnitude of the impact of individual layers of heterogeneity (defined as: ecosystem components that may exert a spatio-temporally uneven influence over the soil surface thermal regime at the decimetre scale) and their interactions on system thermal complexity is largely unknown. This is likely due to limitations in temperature measurement at the scale of interest (decimetre) until recently (Leonard et al 2018). The extent to which combinations of layers induce hotspots or hot moments of varying intensity and longevity, that may induce thermally complex ecosystems that can support biological and biogeochemical diversity or exceed localised tipping points, is not known. It is further unclear how thermal heterogeneity develops from, or can be maximised within, structurally homogenous peatland landscapes, such as within disturbed peatlands, e.g. after fire or harvested peatlands (Turetsky et al 2002).
Due to the number of structural layers and, therefore, the number of ecosystem permutations required, measurement of the impact of each structure and its interaction at a high spatio-temporal resolution can only be determined through numerical modelling. Structural components can be isolated to quantify their individual impact and the impact of their interactions. We simulate the spatio-temporal thermal regime of peatlands using a peatland thermal model that offers the opportunity not only to predict mean peat-surface temperatures but the thermal complexity in space and time (BETA+; Kettridge et al 2013).
The thermal model offers the potential to address our core aim, that is, to determine the influence both of individual layers of heterogeneity (subsurface icelayers, ground-cover vegetation, microtopography, tree, and sub-canopy vascular cover) on surface temperature distributions, and how different combinations of these layers influence peat-surface temperature patterns. We test the conceptual understanding proposed by Leonard et al (2018) that individual layers of system heterogeneity, alone, create sustained heterogeneity in the peat thermal regime, producing sustained hotspots at the peat surface. However, when multiple layers of heterogeneity are stacked, heterogeneity in the thermal regime is reduced because the combined effect of layers results in short-lived low intensity events. To achieve this aim the model is first parameterised and evaluated utilising high spatio-temporal peat temperature measurements of Leonard et al (2018). This novel evaluation using such a high spatio-temporal data sets also identifies model strengths and weaknesses, providing focus to future model development to further our ability to predict the resilience of such heterogeneous and globally significant systems.

Study site
Empirical data was collected from a poor fen (peat-depth ⩾ 3 m) in central Alberta,Canada (55.81 • N,115.11 • W). The site has poor tree cover of Picea mariana with a basal area of 11 m 2 ha −1 , and an average height of 2.3 m (Kettridge et al 2012), which is characteristic of boreal plains peatlands (Wieder et al 2009). The vascular subcanopy includes Rhododendron groenlandicum, Rubus chamaemorus, Chamaedaphne calyculata, Maianthemum trifolia, Oxycoccus microcarpus, Vaccinium vitus-idea and Eriophorum spp. Ground-cover consists of a mosaic of S. fuscum (43%), P. schreberi (9%), Cladina sp., Cladonia sp. (11%) and bare peat (26%) (see figure SI-1 (available online at stacks.iop.org/ERL/16/024002/mmedia)), the combination of which gives the peatland its characteristic, patterned (Sjörs 1961, Foster et al 1983, hummock-hollow features that vary in height by up to 0.4 m within the measurement area, which is consistent with other accounts (Lewis and Dowding 1926). Typically peat freezes each winter, and ice remains until late June or July, or later depending on water levels and winter weather (Petrone et al 2008). Mean ice depth during the study period was 209 mm.

Soil temperature data collection
Soil-surface temperature data was collected using a Silixa Ltd XT fiber-optic distributed temperature sensing (FO-DTS) system buried at 0.02 m depth in the formation of 11 rows spaced 1 m apart in a 10 × 10 m plot. The cable was buried 5 d in advance of data collection to allow for system settling. Temperatures were measured at 0.25 m spatial (along cable length) and 1 min temporal intervals between 06:30 and 20:30 for each measurement day. Temperatures were recorded between 21 May and 3 June 2015 and consist of 4 d of undisturbed conditions, 1 d of trees removed conditions (felled) and 1 d of all vascular vegetation removed (cleared) conditions. Two 10 m lengths of cable were buried outside the plot and remained undisturbed throughout the measurement period as a control. All treatments were applied by hand, using boardwalk between measurement rows, ensuring minimal disturbance to measurement locations. For additional information and photos of the experimental set-up and data-collection methods the reader is directed to (Leonard et al 2018) and associated supplemental material. Standard meteorological data were collected from a 10 m tower approx. 350 m from the plot throughout the measurement period. Weather conditions during the experiment periods were characterized by hot, dry, largely cloudfree conditions, with maximum air temperatures ranging from 25 • C to 28 • C (tables S1-5). No rain fell during the experimental period, except on 31 May and 1 June, between the felled and cleared treatment periods, a small rain event of 13 mm was recorded. Rain-free periods and high surface temperatures subsequently resumed prior to the cleared temperature measurements. At the start of the measurement period, soil moisture (measured using a ML3 Theta probe), the dominant ground cover type (Cladina sp., P schreberi, S fuscum, Bare ground; visual identification) and ice depth (measured using a steel rod) at each measurement and simulation location (0.25 intervals along the length of the FO-DTS cable) were recorded. At the end of the measurement period the tree locations, and photos of crown heights and shapes were recorded (further information may be found in S1-1).

BETA+ model design, parameterisation and evaluation
BETA+ (see S1 for further details) is a finite difference Fourier based model that simulates the 1D thermal behaviour through numerous vertical peat profiles discretised at 0.01 m intervals. The thermal properties of each node are determined based on the peat bulk density and volumetric water content. Volumetric water contents were calculated from the van Genuchten equation, with pore water retention determined from the measured water table depth assuming hydrostatic equilibrium through the peat profile. The volumetric heat capacity was calculated by summing soil constituents multiplied by their respective heat capacities. The impact of ice melt was represented by an increased volumetric heat capacity during periods of melt (see Mckenzie et al 2007a). The thermal conductivity was calculated according to Farouki (1986). Horizontal energy exchange was excluded within the model simulations. Each 1D soil profile was simulated independently of the other. Whilst such horizontal energy exchanges can influence soil temperatures deeper within the soil profile, notably within hummock centres (Kettridge and Baird 2010), the impact on this energy exchanges on near surface peat temperatures simulated here are considered small and were therefore excluded from model simulations for simplicity. Water table depth and soil temperature were measured by a pressure transducer throughout the simulation period and applied to set the water table position and the lower boundary temperature of the model. The water table depth varies amongst soil profiles (by up to 0.41 m) due to the surface micro topography of the simulation plot. Initial peat near-surface temperatures were defined from measured temperatures across the simulation areas at 06:30, at a depth of 0.02 m. Peat profile initial conditions were set using a cubic spline interpolation from measured soil temperatures at the beginning of the simulation period.
The surface boundary of BETA+ was driven by standard meteorological data used to calculate the net short and longwave radiation, sensible and latent heat fluxes. Short and long wave radiation were calculated using a spatially explicit 3D radiation model. Trees within the measurement plot were photographed, and their structures represented as collections of voxels (SI-1). Trees were placed within the 3D landscape and the direct and diffuse shortwave radiation (through time) and sky view factor (constant) at the surface of each peat profile were determined following an approach similar to that presented by Essery et al (2008). The net longwave radiation was calculated in accordance with Stefan-Boltzmann and Beer's law (Aber and Melillo 1991) applying the predetermined sky view factor. Additionally, the model accounts for the influence of micro-topography by incorporating the effect of measured variations in elevation, slope, aspect and shading of the peat surface on the short and long wave radiation balance (Kettridge and Baird 2010). Ice depth measurements were used to inform the spatially varying depth of increased volumetric heat capacity during periods of melt. Latent heat was calculated according to the Penman-Monteith model and sensible heat was calculated according to Newton's law of cooling.
Model performance was evaluated by numerically replicating the experimental setup of the soil temperature data collection, with vegetation layers removed in the same sequence and at the same time (trees removed (felled) then all other vascular vegetation removed (cleared)). The evaluation simulations replicate the experimental design and order. These results are treated separately from the scenarios modelled below.

Model simulation design
Eighty scenarios were simulated, using meteorological data between 06:30 and 20:30 on 24 May (a predominantly cloud free day with air temperature ranging from 4.7 • C to 28.2 • C). Scenarios comprised all permutations of the following options: trees (on/off), shrubs (on/off), microtopography (on/off), surfacecover (on/off), ice cover (on/off), initial surface temperature (on/off). Trees 'on' use spatio-temporally varying solar radiation inputs at the ground surface, as calculated from tree shape, size and location and solar position (described in S1-1). When trees were 'off ' , the mean value for solar radiation was used. When shrubs were 'on' , LAI was assigned from a measured distribution of LAI values obtained at the study site. When shrubs were off, the mean LAI was used. When microtopography was set to off, the surface was represented as one flat surface set at the mean height above water table. When microtopography was set to 'on' , it was represented as a spatially variable surface as measured. Surface-cover was simulated using either a median surface resistance value (median r s , see S1-1) or a surface resistance value selected from a log normal distribution (log-normal distribution r s ) of surface resistance values obtained for each ground-cover type (SI-1). The median surface resistance values allow for a comparison of the influence of heterogeneity associated with differences in surface resistance observed between different surface covers to be assessed. Selecting surface resistance from a log distribution assesses the additional heterogeneity of within ground cover variability in surface resistance. When surface-cover was set to 'off ' , uniform ground-cover was simulated for each of the following conditions; S fuscum (median r s ), S fuscum (log normal distribution r s ), P schreberi (median r s ), P schreberi (log normal distribution r s ), bare ground (median r s ), bare ground (log normal distribution r s ), Cladina sp. (Median r s ), Cladina sp. (log normal distribution r s ). When surface-cover was set to 'On' , the ground-cover was simulated with the species that were present in each measurement location, a mosaic of S fuscum (43%), P schreberi (9%), Cladina sp., (11%) and bare peat (26%) for a median and log-normal r s . Ice cover was simulated as the average observed depth to ice (off) and observed depth to ice at each measurement location (on). Initial surface temperatures were set by either observed surface temperatures (on) or average of the observed surface temperatures (off). Data was simulated for a depth of 0.02 m below the peat surface. The model outputs the temperature for each spatial location at 10 min intervals.

Analysis
Due to the dynamic non-uniform nature of peatsurface temperatures, we used empirical orthogonal functions (EOFs) to provide a compact representation of the intricate and complex data (Rajic 2002) by spatially representing the variance of a dataset. EOF analysis was undertaken in R, using the function 'eof ' from R package 'Spacetime' , version 1.2-2 (Pebesma 2012). Each EOF represents a proportion of the dataset's variance with the first EOF representing the largest proportion of the total variance. This means that the first few EOFs provide insight into most of the variability within a dataset (Dawson 2016). Spearman rank correlation analysis (rho) was applied to compare spatial patterning of the EOFs. A single daily hotspot value (∆T • C) was calculated for each scenario by subtracting the plot mean temperature (figure SI-3) from the greatest mean temperature value at a single location (figure SI-4).

Model evaluation
Change in spatial patterning of modelled near-surface temperatures between treatments is a similar, but subdued representation of measured temperatures ( figure 1 and table 1). The variance explained by the components of the EOF (1-3) are consistent between both modelled and measured data in the undisturbed plot, and effectively represent the observed shifts in response to disturbances. The anomalies demonstrate a consistent and more uniform response in EOF1 compared to EOF2 and 3 for all treatments. EOF2 and 3 demonstrate more polarised responses, at similar spatial scales. The greatest decrease in rho (which indicates the greatest shift in spatial patterning) is observed both in modelled and measured data sets when undisturbed and felled conditions are compared (average decrease in rho value of 0.35 and 0.23 for measured and modelled respectively). A further, much smaller, decrease is observed in both modelled and measured temperatures when felled and cleared scenarios are compared (average decrease in rho values is 0.09 and 0.07 for measured and modelled respectively).
Model simulations regularly overestimate peat surface temperatures under all treatment conditions between 06:30 and 08:30 (figure SI-2). Under undisturbed conditions the model underestimates surface temperatures between 09:30 and 17:30. All model scenarios show best alignment with measured data between 17:30 and 19:30. Cleared conditions show better replication of the median and quartile data, although the first and fourth quartile data underestimate the observed range. High extreme temperatures are best replicated at the start and end of each measurement period but are conservative during the daytime peak in median temperatures (figure SI-2) under felled and cleared conditions.

Increasing the layers of heterogeneity increases the hotspot intensity up to five layers of system heterogeneity
Spatially variable tree, shrub, S. fuscum and mixed ground-cover layers significantly decreased the simulated median temperatures, but spatially variable microtopography, ice layer, species mixture, initial surface temperatures and surface resistances had no significant effect (figures 2, 3 and tables SI-7, SI-8). All layers apart from ice and species mixture had a significant impact on hotspots, with the presence of trees decreasing hotspot extremes, while spatially variable representation of shrubs, microtopography, groundcover species and initial surface temperatures increased hotspot extremes. Surface resistance values drawn from a log-normal distribution increased hotspot extremes relative to simulations using median surface resistance values. All spatially variable layers apart from shrubs, ice and initial surface temperatures showed significant differences in spatio-temporal distribution of surface temperatures.
Peat surface hotspots are zero when the system is homogeneous, i.e. simulations without tree cover, no shrub cover, no microtopography, with no variation in ground-cover species (single species cover), uniform depth to ice, uniform initial surface temperatures and fixed surface resistance values (figure 3). Increasing the layers of heterogeneity increases the hotspot intensity up to five layers of system heterogeneity. The most extreme hotspot intensity is observed ( figure 3) with all the following present: shrubs (on), microtopography (on), S. fuscum only ground-cover, spatially explicit ice depths and starting surface temperatures, and a log-normal distribution of resistance. On average, as the total number of layers of heterogeneity included increases above five, the hotspot intensity decreases (figure 3). The significant changes in simulated hotspots may be due to either changes in the mean temperatures ( figure SI-4), as a result of changes in shading and incidents of direct radiation reaching the ground surface or changes in maximum temperatures (figure SI-3) as a result of variations in slope and aspect, and/or both.

Model evaluation
BETA+ is representing the spatio-temporal peatsurface temperature dynamics, replicating observed changes in spatial patterning of EOF1 between undisturbed, felled and cleared conditions (figure 1, table 1). Although the simulated EOF extremes are lower than measured, this likely reflects the conservative extreme values simulated by the model compared to the measured data (figures 1 and SI-2). Further, the decrease in average rho values between undisturbed and felled is also well replicated by the model when cleared. However, greater variation in spatial patterning (rho (r)) of EOF1 is observed within simulated than observed undisturbed temperatures (range of 1-0.74 compared to 1-0.94; table 1). This may result from variations in the cable burial depth (±1.5 cm). Cable buried deeper will be consistently cooler and cable buried shallower consistently warmer during the simulated daytime growing season temperatures. Any such bias in the measured peat temperatures will impact the ranking of values in the spearman rank correlation analysis, providing stronger correlations between the patterns in EOF1 derived from measured undisturbed data (table 1)   . Despite these limitations, and the evaluation of the model over a period of a few days due to the spatial intensive nature of the evaluation data set, the model replicates changes in spatial variance with treatment similar in extent to the observed peatland, allowing confident determination that the model replicates observed shifts in peatland thermal dynamics. Therefore, this model for the first time, allows confident separation of individual system structural components on highly variable spatio-temporal near surface soil temperatures.

Influence of heterogeneity on surface thermal regime and system functioning
Simulations broadly support the hypothesis of Leonard et al (2018) that peat surface hotspot intensity is linked to ecosystem heterogeneity (figure 2). However, maximum hotspots are observed when five layers of heterogeneity are present (figure 2) rather than one as previously hypothesised. The 9.1 • C range in hotspot temperature when five layers are present suggests that not all system layers of heterogeneity have an equal influence. The hypothesis that the maximum hotspots are observed only when one layer of heterogeneity is present (Leonard et al 2018) is likely true if all layers have equal influence over the thermal regime and the heterogeneity of each layer is independent from the next. However, co-dependence in the observed system layers is likely (e.g. between trees and sub-canopy vascular species as they compete for resources such as light) resulting in layers of heterogeneity with varying influences over system hotspots. Ice depth variation and groundcover type are the only layers that do not significantly influence peat surface hotspots (figure 2 and tables SI-7, SI-8). As a result, the increased heterogeneity of the system magnifies the intensity of hotspots through the additive effects of the different elements over one another. A reduction in hotspot intensity where layers counter each other's effects is only evident in very complexity multilayer systems, and thus likely in mature fully diverse peatland ecosystem. This should be consisted in the management and restoration of peatlands. In restoration of degraded systems or formation of new systems, building different layers of heterogeneity support the thermal diversity and potential wider ecosystem diversity. But also, disturbance of mature systems that removes such complexity will likely increase uniformity of the thermal systems and potentially of wider associated ecosystem characteristics.
Peatland system functioning is a balance of spatially varying process feedbacks that include; hydrological processes, ecological succession and competition, productivity and decomposition (Waddington et al 2015) Critically, results show that changes in the spatio-temporal dynamics may occur, potentially inducing an overall shift in system functioning, without a significant change in median temperature (figure 2 and tables SI-7, SI-8). Notably, this is observed to occur in response to the addition of microtopography, the inclusion of mixed Figure 2. Each of the three tiles contain data relating to the daily median temperature ( • C), daily hotspots (∆T • C) and changes in spatial distributions (rho (r) values closest to one show least change from the undisturbed system). Each component on the x-axis represents data from all scenarios, with the different colours representing the presence/absence/type of layer heterogeneity tested. (a) Trees (data from all scenarios with trees on vs data from all scenarios with trees off), (b) shrubs (data from all scenarios with shrubs on vs data from all scenarios with shrubs off), (c) micro-topography (data from all scenarios with micro-topography on vs data from all scenarios with microtopography off), (d) mixed species (data from all scenarios with mixed ground-cover on vs data from all scenarios with single species ground-cover), (e) species (data from all scenarios comparing each ground-cover type), (f) ice (data from all scenarios with ice on vs data from all scenarios with ice off) (g) initial surface temperature (observed surface temperatures (on) vs average of the observed surface temperatures (off) and (h) resistance (data from all scenarios with log-normal distribution of surface resistances (rs) vs data from all scenarios with fixed value surface resistances (rs)). Boxplots represent the median, 25th and 75th percentile (whiskers: smallest and largest observed value that's less than or equal to the lower/upper hinge ±1.5 × IQR). Asterisks denote significant differences between each layer option, calculated using Wilcoxon rank sum test ( * = <0.05, * * = <0.005, * * * = <0.0005; table SI-6).
species and as a result of within species variability in surface resistances. Differences in the drivers of such peat-surface thermal patterning and dynamics occur between different peatland systems but may also vary in time as result from a range of disturbances within a peatland that impact the layers of heterogeneity discussed. These may include anthropogenic disturbances (such as felling, exploitation for energy etc), biological disturbances (e.g. insect infestations or plant disease) and climate driven disturbances (increased fire frequency, shifts in species distributions and composition, water availability, extreme weather events etc). Such disturbances will shift rates and locations of thermally driven processes such as productivity, species competition, decomposition, and evapotranspiration. As a result, the system is placed under a heterogeneous stress because hotspots have shifted locations and may have dis-proportionate effects on system process rates (Johnstone et al 2016). Locations where inherent resilience to increased peat surface temperatures had developed, e.g. microbial acclimation to elevated temperatures (Bradford et al 2008, Kaiser et al 2014, no longer align with the new peatsurface temperature signature. If spatial changes result in process rates that breech thresholds and tipping points, they could irreversibly shift the balance of key system feedbacks (Rietkerk and van de Koppel 2008, Belyea 2009, Waddington et al 2015 and cause major shifts in system functioning (Johnstone

Conclusion
Our results clearly show changes in spatial distributions of temperatures both with and without changes in median values. Models used to predict changes to system functioning rarely consider the heterogeneous nature or interactions of disturbances on systems where significant changes in the average values of process drivers are observed, let alone heterogeneous shifts where no significant changes in average process drivers are observed. To fully understand system functioning and its response to disturbance, a more comprehensive understanding of spatio-temporal responses to heterogeneous stresses is recommended. Researchers should consider the spatio-temporally heterogeneous nature of any disturbance and the likely impact of this on the spatio-temporal dynamics of balanced system feedback mechanisms.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).
Worrall and three anonymous reviewers for their feedback on this work.