Modelling and upscaling ecosystem respiration using thermal cameras and UAVs: Application to a peatland during and after a hot drought

Abstract Field-based thermal infrared cameras provide surface temperature information at very high spatial and temporal resolution and could complement existing phenological camera and spectral sensor networks. Since temperature is one of the main drivers of ecosystem respiration (ER), field-based thermal cameras offer a new opportunity to model and upscale ER in unprecedented detail. We present such an approach based on manual chamber CO2 flux measurements and thermal imagery from a tower-based camera and from Unmanned Aerial Vehicle (UAV) flights. Data were collected over two growing seasons, including the hot drought of 2018, for the two main vegetation microforms (hummock and hollow) of a hemi-boreal peatland in Sweden. Thermal imagery proved suitable for modelling ER in this ecosystem: ER model accuracies were similar when air, soil or surface temperature measurements were used as input. Our findings allowed us to upscale ER using UAV-derived thermal images and we present maps of ER at sub-decimeter resolution (


Introduction
Respiration is the production of carbon dioxide (CO 2 ) by living organisms as they grow, live and decompose and ecosystem respiration (ER) is the sum of respiration from all vegetation and soil fauna within a specific ecosystem. Monitoring ER is essential because it represents a CO 2 source to the atmosphere 10 times larger than anthropogenic CO 2 emissions and has a direct impact on the amount of carbon (C) that remains stored in the biosphere (IPCC, 2014;Yuan et al., 2011). Temperature is one of the main factors affecting ER rates (Atkin and Tjoelker, 2003). Therefore ER is often simply modelled based on an exponential or polynomial fit with temperature, although soil moisture, productivity and substrate availability also affect ER (Davidson and Janssens, 2006;Kirschbaum, 1995;Lloyd and Taylor, 1994).
Since ER is the sum of multiple above-and below-ground processes, there are ongoing discussions of what temperature metric (air, soil or surface temperature) is most appropriate to use in ER modelling (Lasslop et al., 2012;Wohlfahrt and Galvagno, 2017). One of the major limitations of temperature measurements is their coarse spatial resolution, with air and soil temperature sensors usually limited to one or a handful of locations around a study site. Recent studies have highlighted this limitation, by showing that microclimate can vary widely within an ecosystem, thus leading to large deviations from standardized air and soil temperature measurements at a single location (Lembrechts and Lenoir, 2020;Zellweger et al., 2019). For example, leaf surface temperatures may differ by >40 • C from air temperature (Still et al., 2019) whilst maximum air temperature can vary by 20 • C over a few hectares (Maclean, 2020).
Field-based thermal infrared cameras offer a new approach for modelling and upscaling ER because they capture images of surface temperature (where 'surface' represents the top layer of the vegetation canopy or bare soil) at high spatial and temporal resolution. A thermal image can contain >10 5 pixels, each of which provides a separate temperature measurement. Satellite-derived land surface temperature (LST) data provide global coverage but due to their coarse resolution (pixel size: 60 m to 1 km, temporal resolution: 1 to 16 days) they cannot be used to study microclimatic variability. Ground, Unmanned Aerial Vehical (UAV) and airborne thermal cameras however can measure temperature at an extremely high spatial resolution (pixel size: <10 cm to ~50 m, temporal resolution: minutes to months), which is similar to, or finer than, the spatial scale at which measurements of ecosystem carbon fluxes are conducted. These technologies can thus be used to investigate how variations in microclimate across an ecosystem affect carbon fluxes and other biological processes (Still et al., 2019). Although satellite and airborne LSTs have previously been used to model ER and soil respiration, including in peatlands, (e.g. Anderson et al., 2008;Jägermeyr et al., 2014;Schubert et al., 2010), the use of images from thermal cameras deployed on towers or UAVs to model and upscale ER has, to our knowledge, not been tested in any ecosystem. Mapping the spatial variation of ER across an ecosystem is valuable because such maps could improve our understanding of flux tower measurements and their potential biases by showing how ER varies within the flux tower footprint (i.e. source area) and beyond.
Peatlands are a prime testing-ground for evaluating whether there are benefits to using thermal imagery (i.e. high spatial resolution temperature data) to model ER. Peatlands are not only globally significant carbon stores (Yu et al., 2010), but in many cases their unique microtopography creates a spatially heterogeneous ecosystem (Waddington and Roulet, 1996). The characteristic hollow and hummock landscape of a boreal peatland creates variations in hydrology, microclimate and vegetation species composition that can lead to large differences in carbon fluxes within a peatland (e.g. Cresto Aleina et al., 2015;Dieleman et al., 2017;Sullivan et al., 2008). Furthermore, peatland surface temperature exhibits significant spatial variation. Mölder and Kellner (2002) reported surface temperature fluctuations of >10 • C along a 100 m transect in a Swedish mire, while Leonard et al. (2018) observed temperature ranges exceeding 25 • C within a 100 m 2 area of a Canadian fen. Such large spatial temperature variation would have a significant impact on modelled ER (ER mod ) estimates due to the exponential relationships used, whereby a small temperature change produces a large change in ER mod . As a result, current ER estimates which are based on only one or few spatially distributed temperature measurement points may underestimate the variability of ER and fail to capture hotspots of carbon release. Upscaled carbon flux estimates may thus be biased or unrepresentative of the whole ecosystem. This study investigates, at a hemi-boreal peatland, what improvements high resolution surface temperature data from tower-and UAVmounted thermal cameras could bring to ER modelling compared to using sparse air and soil temperature data. Since the thermal camera data is spatially resolved, we start by assessing whether there are spatial differences in temperature and ER, focusing on the two main peatland vegetation communities (hummocks and hollows) at our study site. Secondly, we evaluate the performance of ER mod driven by surface temperature data compared to using point measurements of soil and air temperature. Thirdly, we demonstrate the potential advantages of UAVderived surface temperature data for upscaling ER, by showing how it can be used to create high resolution maps of ER across our study area, whilst also noting the limitations of such an upscaling approach. Since our fieldwork was conducted over two growing seasons, including an extreme drought year, we further assess how the drought affected surface temperatures and ER.

Methods
The methods of this study include four main steps: (1) collection of observational data (temperature, ER and ancillary data), (2) ER modelling based on the observed ER and temperature data, (3) analysis of the observational data and modelling results, and finally, (4) ER upscaling using data collected using UAVs and the ER models developed in step (2).

Study site
Data collection took place at Mycklemossen, a hemi-boreal fen with oligotrophic bog-like vegetation in southern Sweden (58 • 21'N, 12 • 10'E, Fig. S1) in 2018 and 2019. The site is part of the Skogaryd Research Catchment and Swedish Infrastructure for Ecosystem Science network. The long-term (1981-2010) mean annual air temperature and total precipitation were 6.9 • C and 901 mm, respectively, at the closest national monitoring stations. Due to the strong precipitation gradient between the site and the coastline, we averaged precipitation data from two national monitoring stations (Vänersborg, 10 km east of and 30 m lower elevation to Mycklemossen, and Uddevalla, 14 km west of and at similar elevation to Mycklemossen; SMHI, 2019). Temperature data was only available from Vänersborg. We present precipitation data for 2018 and 2019 from the two national monitoring stations, and air temperature data that was recorded at Mycklemossen (see below for details). During the late spring and summer of 2018, the site experienced a 'hot drought', i.e. hot and dry conditions. To assess the severity of the drought in 2018, we present values of the Standardized Precipitation Evapotranspiration Index (SPEI), a drought index based on precipitation and potential evapotranspiration with a mean of 0 and standard deviation of 1 (Vicente-Serrano et al., 2010). The SPEI was retrieved at a 1 • spatial resolution for the site in 2018 and 2019 from the Global Drought Monitor Project (Vicente-Serrano et al., 2020).
The site is composed of a mosaic of drier hummocks and wetter hollows, with local peat deposits extending down to 6 m below the surface (Wallin et al., 2015). The hummocks are dominated by Eriophorum vaginatum, Calluna vulgaris and Erica tetralix whereas the hollows consist of different Sphagnum species, mainly S. rubellum, S. fallax and S. austinii as well as Rhyncospora alba.

Ecosystem respiration
We conducted manual ER measurements using a non-steady state through-flow chamber system (Livingston and Hutchinson, 1995). Six stainless steel collars (7 cm depth, 47 cm x 47 cm) were installed in May 2018, divided evenly between hollow and hummock areas. The collars were placed in the field of view of a thermal camera mounted on a mast ( Fig. 1 and Section 2.1.3). We used six collars because this was the maximum number that could be accommodated within the centre of the camera's field of view (to avoid vignetting effects, see Section 2.1.3) without disturbing each other and so as to minimize trampling during sampling. The area seen by the thermal camera is representative of the vegetation and conditions across the whole fen (Fig. 2). The collars were classified as hollow or hummock when the vegetation was dominated by Sphagnum sp. and Rynchospora alba or by vascular vegetation, respectively (Fig. 1). Since our collars include vegetation, we use the term 'ecosystem respiration' (ER; as opposed to only 'soil respiration') to refer to our CO 2 chamber measurements and subsequent modelling.
ER was derived manually with an opaque chamber (47 cm x 47 cm x 30 cm) covered in aluminium foil. The CO 2 concentration was measured at 1 Hz with an Ultraportable Greenhouse Gas Analyzer (UGGA; Los Gatos Research, San Jose, USA). A small fan circulated air inside the chamber and an external pump was used to aid circulation between the chamber and the UGGA (located 30 m apart). The chamber closure time was 5 min. During chamber closure, the standard deviation (σ) of the air temperature (as measured by the UGGA) was on average 0.6 • C. The height of the vegetation surface relative to the height of the collar was accounted for when calculating the chamber volume. Monthly campaigns from June-September resulted in 516 ER measurements (324 in 2018, 192 in 2019). To derive CO 2 fluxes from the concentration measurements, the first 30 seconds of each measurement were discarded to avoid disturbances caused by the chamber being placed on the collar (Davidson et al., 2002). Next, five linear regressions (CO 2 concentration against time) were fitted to the data over 3 min windows, with an interval of 10, 20 or 30 sec between the start of each window. The slope of the regression with the highest R 2 above 0.9, as well as half-hourly mean air temperature and pressure measurements from the site, were used to calculate the flux using the ideal gas law. ER measurements where all regressions had R 2 <0.9 (3% of the measurements, either cold periods with low fluxes or when the chamber was not completely sealed) or where corresponding abiotic data were unavailable (9%), were discarded. All other measurements within the same round (one 'round' is one measurement of each of the six collars) were also discarded to ensure an equal distribution of hollow and hummock measurements. As a result, 414 measurements were used in the analysis (228 from 2018 and 186 from 2019).

Temperature
We used three types of temperature measurement: air, soil and surface. Air temperature and humidity were measured at 1 Hz with an HC2S3 sensor (±0.1 • C accuracy, Campbell Scientific, Logan, UT, USA) 65 m away from the collars at 2 m height in a ventilated, radiationprotected housing. Mean air temperature was calculated using all measurements in the 15 minutes before and after each flux measurement to ensure consistency with the surface temperature measurements (see below). Soil temperature was measured with a handheld digital soil thermometer (±1 • C accuracy, model 6310, Spectrum Technologies, Inc., Aurora, IL, USA) 6 cm below the surface, once after each ER measurement, outside the collar (within 10 cm of the collar). The mean range of soil temperature during a measurement round was 2.3 • C and 1.1 • C for the hummocks and hollows, respectively.
Surface temperature was measured using a FLIR A65 thermal infrared camera (13 mm lens, 640 × 512 pixel resolution and 7.5-13 µm spectral range; FLIR Systems, Inc., Wilsonville, USA). The camera accuracy is ±5 • C or 5% of the reading and the thermal sensitivity is <0.05 • C at 30 • C. The camera was installed on a tower 9 m above the surface of the fen, facing nadir (referred to as the 'fixed' camera). A calibration plate (see Kelly et al., 2019 for details) and the collars were placed in the centre of the camera's field of view (FOV) to minimize the impact of vignetting (reduction in temperature values from the centre to the edge  of the camera's FOV; Mitchell, 2010; see Fig. 1). Initially images were captured every 5 min, but the frequency of image capture was increased to every minute during ER measurements after the first half of the 2018 field season. The increased image frequency had no significant impact on the variability of the average surface temperature (mean σ = 2.23 • C and 2.12 • C, before and after the increase in image frequency, respectively).
Thermal image data processing followed the method described in Aubrecht et al. (2016) whereby 14-bit raw data were converted to surface temperature measurements using air temperature, sky temperature, humidity, distance to target and emissivity input data. In addition, a flatfield correction was performed using a Gaussian filter to minimize the impact of vignetting (Leong et al., 2003). The images were also corrected for temperature drift (noise related to changes in camera temperature) by performing a one point offset correction using the known surface temperature of the calibration plate (FLIR, 2016). Gaps in the calibration plate surface temperature (due to thermistor failure) were filled using a linear regression based on the surface temperature of a similar, nearby calibration plate (model RMSE = 2.13 • C). The reduction in transmission of infrared radiation by the window of the camera housing (transmittance = 88%) was also accounted for. As a result, the mean difference between the surface temperature of the calibration plate measured by a thermistor (Model 107, ±0.4 • C accuracy, Campbell Scientific Inc., Logan, USA) and that measured by the camera was 1.13 • C (σ = 0.80 • C) during the growing season. Vegetation emissivity was fixed at 0.97, an average of literature values for peatland vegetation (Goïta et al., 1997;Kettridge and Baird, 2008;Mölder and Kellner, 2002). After processing, mean surface temperatures were extracted from the pixels inside each collar for images within a 15 min period before and after each flux measurement. Each collar represented ~860 pixels.

Water table depth
Water table depth was recorded at 0.1 Hz using CS450 pressure sensors (Campbell Scientific, Logan, UT, USA) at two points (one hummock and one hollow) within 60 m of the collars. Since the area is predominantly flat, these sensors provide a sufficient estimate of the temporal fluctuations in the water table at our collars.

Modelling of ecosystem respiration
We modelled ER (ER mod ) using the Global Polynomial Model (GPM) from Heskel et al. (2016;Eq. (1)) using surface, air or soil temperature as inputs (ER mod_surf , ER mod_air , and ER mod_soil , respectively). The GPM is used to model ln(ER mod ) but for clarity we present the results transformed to ER mod . We use the GPM because it successfully captures the decline in respiration rate at high temperatures (see Heskel et al., 2016), a biological phenomenon (Atkin and Tjoelker, 2003). Moreover, it was developed using surface temperature, which we also use here.
where T = temperature and a, b and c are free parameters. Heskel et al. (2016) determined that b = 0.1012 and c = 0.0005 based on their global dataset of leaf respiration and temperature and suggested that values of a should be parameterized based on the Plant Functional Type of interest. We tested the GPM model with one (a; 'GPM1') or three (a, b, c; 'GPM3') free parameters. We fit the models using Leave-One-Out-Cross-Validation (LOOCV), calculating the mean of the model coefficients from each CV model fit for the final model. We used Akaike's Information Criterion (AIC) to decide whether the improved fit provided by an increased number of free parameters outweighed the increase in model complexity. The model with the lowest AIC is regarded as the most parsimonious model. In this case, the AIC of GPM3 was similar or lower across both years and all temperature variables than GPM1 (Table S1). Therefore, we continued our analysis using only GPM3, which is hereafter referred to as 'GPM'. We present the GPM fit with 18 different parameter sets corresponding to separate temperature-ER relationships for 3 vegetation categories (hummock, hollow and both combined), each driven by 3 temperature metrics (soil, surface and air) for two years of observations (2018 and 2019).

Analysis of observed and modelled data
We began by testing for significant differences in mean surface, soil and air temperature as well as ER between hummocks and hollows (for the times when ER was measured). We used t-tests or Wilcoxon tests, depending on whether samples were normally distributed according to a preliminary Kolmogorov-Smirnov test (Kaltenbach, 2012). We also tested for significant differences in these variables aggregated on a monthly basis. In addition, we examined the relationship between water table depth, ER and all three temperature metrics, and conducted partial least squares regressions to understand how the covariation of water table depth and temperature impacted ER.
To assess whether the high spatial resolution surface temperature data provided any improvement in ER mod accuracy over air or soil temperature, the models were validated using LOOCV. Due to the hot drought in 2018, we modelled ER in 2018 and 2019 separately. To assess the impact of vegetation heterogeneity on ER mod , for both years, we investigated whether the observed ER fluxes were best captured using a single model or using separate models for hummocks and hollows by comparing the model normalized root mean square error (NRMSE, normalized using the mean of the observed ER) and σ.
To assess the impact of the choice of temperature metric (surface or air) on ER mod , we compared the mean growing season ER mod_surf, and ER mod_air (based on all available air or surface temperature data between June-September; soil temperature data was only available for times when flux measurements were performed). We estimated mean growing season ER mod , rather than the sum, as technical issues with the thermal camera led to gaps in the surface temperature data that could not be gapfilled. In 2018, the surface temperature data were missing for parts of July and August while in 2019, data are mostly missing in June. To account for the impact of the uncertainty in the temperature measurements on ER mod , we also calculated mean growing season ER mod using temperature data with an added random error that was proportional to the manufacturer stated error: ±5 • C and ± 0.1 • C for surface and air temperature, respectively. The error was randomly sampled from a normal distribution (mean = 0, σ = manufacturer error) and added to the original temperature values. Finally, the temperature and ER mod fluxes for each hour of the day were averaged and plotted to show the diurnal cycle of temperature, ER mod_surf, and ER mod_air .

Upscaling to ecosystem scale
We upscaled ER mod_surf over a larger area of the fen (120 m x 110 m) based on surface temperature data collected with an Unmanned Aerial Vehicle (UAV). We used a Solo (3DR, Berkeley, USA) and an Explorian 8 (Pitchup AB, Mölndal, Sweden) quadcopter. Details of the UAV flight protocol are outlined in Kelly et al. (2019). We conducted 22 flights during the 2018 and 2019 growing seasons to collect surface temperature data. Flight duration was up to 12 minutes with average speeds of 2-6 ms − 1 while flying at 50 m above ground level, yielding a spatial resolution <7 cm. Surface temperature data were collected using a FLIR Vue Pro 640 and a FLIR Vue Pro 640 R with ThermalCapture Calibrator (TeAx Technology GmbH, Wilnsdorf, Germany) in 2018 and 2019, respectively. The cameras have a sensor resolution of 640 × 512 pixels, a 13 mm lens and a spectral range of 7.5-13 μm. The Vue Pro 640 has an approximate accuracy of ±5 • C, provides raw data in 14-bit tiff files and is radiometrically calibrated by the user (cf. Kelly et al., 2019), whilst the Vue Pro 640 R is radiometrically calibrated by the manufacturer to ±5 • C and provides data in radiometric jpeg files (including 14-bit raw data). The ThermalCapture Calibrator improves the accuracy of surface temperature measurements by up to 70% by providing an external shutter with a known temperature that is used to improve the automated, internal camera calibration (TEAX Technology, 2020). In addition, we conducted one UAV flight in June 2017 using an RGB (A6000, SONY, Weybridge, UK) and multispectral camera (RedEdge, Micasense Inc., Seattle, USA) to collect images for a land cover classification of the fen. All images were processed in Metashape Professional (Version 1.4.1, Agisoft LLC, St. Petersburg, Russia) to produce orthomosaics for each flight.
To identify hummocks and hollows (Fig. 2), we classified the RGB and multispectral orthomosaics of the fen using random forest classification (Breiman, 2001) with the TreeBagger function in MATLAB (50 trees and default input parameters; MathWorks, Natick, MA, USA). The following 11 bands were used as input: RGB (all bands), multispectral (red, green, blue, red-edge and near infrared), Normalized Difference Vegetation Index (based on both the red and red-edge bands) and a digital surface model (produced in Metashape Professional).
To train and validate the classification, vegetation species data were collected in the field along 12 transects divided equally between a training and a validation area (30 m x 30 m), on the east side of the fen. These areas where chosen because they were easily accessible from the edge of the fen (minimal trampling) and the vegetation was representative of the whole fen. We used the data collected in the training area (460 pixels) and a visual assessment of the RGB orthomosaic to manually label >30 million pixels (3.5 cm resolution) as a training dataset for the classifier using ENVI image analysis software (Harris Geospatial Solutions, Inc., Broomfield, CO, USA). Labelling the additional pixels improved the classification, particularly the representation of the tree class. The classification was validated using the data (580 pixels) from the validation area (total accuracy = 81%). After applying a 5 × 5 modal filter to reduce salt-and-pepper noise, the classification yielded a total accuracy of 84% when validated using the fieldwork data (see Table S2 for full confusion matrix). We used a pixel-based classification approach (instead of an object-based approach) because the vegetation communities in the fen did not produce distinctive objects and it would have required ground-truth training and validation data to be collected at the object-level which would have been much more labour intensive.
We present thermal orthomosaics from five UAV flights (see Table 1). The temperatures from each orthomosaic were used as input into the ER model fits driven by surface temperature for hollows and hummocks (see Section 2.2). The locations of the hollows and hummocks were identified using the land cover classification. As a result, we could map ER mod_surf over the area covered by each thermal orthomosaic. The classification and orthomosaics presented here include the area observed by the fixed thermal camera and the collars used for the ER measurements (Fig. 2). To assess the impact of the uncertainty in the UAV thermal camera temperature measurements (±5 • C), we repeated the same uncertainty analysis as described in Section 2.3, and calculated the mean ER mod_surf for each orthomosaic, both including and excluding the additional random error in the surface temperature data. To analyse the impact of the spatial heterogeneity of ER, we further compared the above mean ER mod of each orthomosaic derived by either accounting for (i.e. ER mod_surf ) or by neglecting (ER mod_surf_all ) the proportions of hummocks and hollows and their seperate GPM model parameter sets.

Climate data
In 2018, large parts of Sweden experienced a hot drought (hot and dry conditions). Monthly mean air temperature at Mycklemossen reached or exceeded the climate normal range between May and November by up to 2.5 • C (Fig. 3a). Monthly precipitation sums were below average in late winter and from May to July, with extremely dry conditions in July when only 20 mm of rain was recorded, compared to the long-term average of 80 mm (Fig. 3b). As a result, the water table was low during most of the growing season, dropping to − 0.28 m and − 0.45 m below the surface for the hollows and hummocks, respectively, at the peak of the hot drought (Fig. 3c).
In 2019, climatic conditions were generally closer to the long-term average but with monthly air temperatures down to 3.5 • C lower than the long-term average in April and May (Fig. 3a). Total monthly precipitation was above the long-term average before and after the growing season, with a short dry-spell in April (Fig. 3b). The water table level declined during the growing season, but for a shorter duration than in 2018, reaching minimum levels of -0.22 m and -0.38 m for the hollows and hummocks, respectively (Fig. 3d).
The 3-monthly SPEI at the site reached a minimum of -2.57 in July 2018, at the peak of the hot drought, and -1.49 in June 2019 (Vicente-Serrano et al., 2020 ). Due to the severity of the hot drought during the peak of the growing season in 2018, we analysed our ER measurements and modelled ER separately, for 2018 and 2019.

Statistical differences between hummock and hollow temperatures and ER
To assess the spatial variability of temperature and ER across the fen, we tested for significant differences in temperature and ER between the hummocks and hollows. We found no significant differences in surface or soil temperature between the hummocks and hollows in 2018 or 2019 (for the times when ER chamber measurements were conducted, Fig. 4a). Monthly means of surface and soil temperature, however, differed significantly between hummocks and hollows during the hot drought in 2018 (Table S3). The mean surface temperature of the hollows was 3 • C warmer than the hummocks, whereas the mean soil temperature of the hollows was ~1 • C cooler than the hummocks in July 2018. During that month, mean surface temperatures reached 41 • C whereas air and soil temperatures were <28 • C. Mean growing season surface, air and soil temperature of both vegetation communities combined were 2.0 • C, 2.8 • C and 1.8 • C warmer in 2018 than in 2019.
ER was significantly different between the two vegetation communities in both 2018 and 2019. In 2018, the hollows had higher mean ER whereas in 2019, the hummocks had higher mean ER (Fig. 4b). Mean ER declined during the hot drought for both vegetation communities, but the hummocks experienced a larger decline (48%) compared to the hollows (15%).

ER mod and choice of driving temperature
We used the observed ER fluxes to fit the GPM (Eq. (1)) to 18 ER mod expressions based on different parameter sets to account for the hollows and hummocks (separately and combined) and three different driving temperatures (surface, air or soil), for both 2018 and 2019 (Fig. 5). The model parameters are listed in Table S4. The choice of temperature metric produced only small differences in model NRMSE, although differences in R 2 were more pronounced (Table 2). In 2018, ER mod_surf produced the best fit whilst in 2019 ER mod_air produced the best fit ( Table 2). The most notable difference between model fits driven by the three temperature metrics was the change in temperature sensitivity of ER mod . Since surface temperature spanned a wider range than air or soil temperatures (40 • C compared to 24 • C and 23 • C, respectively, over both years), ER mod_surf had flatter slopes (i.e. was less temperature-sensitive) compared to ER mod_air or ER mod_soil (Fig. 5). The difference in the temperature sensitivity was especially marked during the hot drought in 2018, when hollow ER mod_surf and ER mod_air had a flattening curve or even a decline in ER mod at high temperatures whereas hollow ER mod_soil showed continuous increases with temperature. When both vegetation communities were combined, ER mod based on all three temperature metrics showed a strong response to the 2018 hot drought: ER mod stopped increasing at high temperatures. Using separate model fits for hollows and hummocks did not consistently improve the model fit compared to using a single model for both vegetation communities combined ( Table 2). The use of separate ER model fits is, however, justified by 1) the significant differences in ER (see Section 3.2) and the model parameters (Table S4) between the two vegetation communities; 2) the separation (i.e. no overlap) of hollow and hummock ER mod over certain temperature ranges (Fig. 5); and 3) the model shape of hollow ER mod being closer to exponential, whereas hummock ER mod showed a faster decline in temperature-sensitivity at high temperatures (Fig. 5).
We also investigated the impact of using air versus surface  temperature data on diurnal ER mod for the whole growing season (Fig. 6). Although ER mod_surf and ER mod_air have similar model accuracy, mean growing season ER mod_surf was between 8-16% lower (depending on year and vegetation community) than ER mod_air (Table S5). The difference is due to surface temperature declining more rapidly than air temperature in the evening plus lower nighttime surface temperatures. In 2019, ER mod_air also increased more rapidly than ER mod_surf in the morning because of the higher temperature sensitivity of ER mod_air compared to ER mod_surf . Uncertainty in the temperature measurements had little impact on mean growing season ER mod (Table S5). The uncertainty (±0.1 • C) in the air temperature measurements made no difference to mean growing season ER mod_air . When the uncertainty in the surface temperature measurements (±5 • C) was accounted for, mean growing season ER mod_surf was <5% higher (depending on year and vegetation community, mean = 3.9%) compared to when it was calculated using the original surface temperature data. Despite the additional error in the input temperature data having a mean of 0 • C, it always caused increases in mean ER mod_surf with the error compared to the original mean ER mod_surf . This increase can be explained by the ER model shapes (Fig. 5) relative to the mean growing season surface temperatures (15 -17 • C): when the input temperature increased due to the additional (1)) separately for hummocks (Hum) and hollows (Hol; a-f) and for both vegetation communities combined (g and h). The GPM is fitted for 2018 (a, c, e and g) and 2019 (b, d, f and h) data using three temperature metrics: surface, air and soil temperature (see Table S3 for model coefficients). Model σ (calculated from all model fits during LOOCV) is too small to be visible. error, it caused larger changes in ER mod_surf than when it decreased.

ER and water table depth
Observed ER increased as the water table dropped, although ER for both vegetation types was lower in 2018 (a warm and dry year) than in 2019 (Fig. 7). ER was more sensitive to changes in the water table in 2019 compared to 2018. Furthermore, water table depth explained more of the variation in hollow ER (adjusted R 2 = 0.78 and 0.69 in 2018 and 2019, respectively), compared to hummock ER (adjusted R 2 = 0.30 and 0.48, Fig. 7).
The relationship between ER and temperature or water table depth (Figs. 5 and 7), may be partially explained by the covariation of surface temperature with water table depth. Declines in the water table were associated with increases in surface temperature (as well as with increases in air and soil temperature, not shown). For both years and vegetation communities, Pearson's r < -0.72 when water table depth was correlated against surface, air or soil temperature. A partial least squares regression of surface temperature and water table depth against ER showed that the first component of the regression explained >70% and >50% of the variance in ER for hollow and hummock ER, respectively, across both years.

Upscaling ER mod to ecosystem scale
UAV-derived surface temperature data offer the possibility to upscale ER mod_surf at very high spatial resolution (pixel size < 7cm). The resulting maps show large variations in ER mod_surf over the study area which are a function of the vegetation community dependent GPM fits, surface temperature variation and the hot drought conditions in 2018. For example, Figure 8a shows that the fen surface temperature varied by more than 15 • C across a 0.5 ha area during the hot drought. The surface temperature was related to the hummock and hollow microtopography, with maximum temperatures occurring on the south sides of hummocks and in some hollows where the top Sphagnum layer completely dried out. The resulting ER mod_surf map (Fig. 8c) shows the lowest ER mod_surf in hummock areas that experienced the highest surface temperatures, whilst the highest ER mod_surf values occur in hollows. In 2019, the surface temperature (Fig. 8b) showed less variability (σ = 2.5 • C compared to σ = 4.2 • C in 2018), which could be due to the higher water table. The  higher temperature variation during the hot drought is supported by the data from the fixed thermal camera, which showed significant differences in surface temperature between hummock and hollows in July 2018 but not in July 2019 (Table S3). The spatial patterns in ER mod_surf in 2019 (Fig. 8d) are driven by both temperature and vegetation community: warmer areas with a higher proportion of hummocks in the southern section of the orthomosaic have higher ER mod_surf than the cooler areas dominated by hollows in the northern section. Fig. 9 shows a larger area of the fen (120 m x 110 m), at different times of the same day in September 2019. In (a-b) and (d-e) the effects of changes in sun azimuth angle on the surface temperature are evident. In (a) and (d), the sun shines from the SE, whereas in (b) and (e), the sun shines from the SW. As a result, the hummock sides facing SE or SW, respectively, are warmer, and thus have higher ER mod_surf than the surrounding areas. In (c) and (f) the sun has set and the surface has cooled rapidly, resulting in low ER mod_surf values. Hence, although the hollow vegetation is warmer than the hummock vegetation, the ER mod_surf is low for both vegetation communities (cf. Fig. 5b).
We assessed the effect of the surface temperature measurement uncertainty (±5 • C) on the mean ER mod_surf of each orthomosaic (Table  S6). Mean ER mod_surf changed between -6% to +7% when the error was included. In all orthomosaics from September 2019 (Fig. 9), mean ER mod_surf increased with the additional uncertainty. The low mean surface temperatures of these orthomosaics (between 2-13 • C) caused the increases in the input temperature data due to the additional uncertainty to have a larger impact on ER mod_surf than decreases in the input temperature. On the other hand, for the orthomosaic in Fig. 8b, there was no difference in ER mod_surf with or without the additional uncertainty because surface temperatures were high (mean = 36.4 • C), and the model slope (Fig. 5b) caused both increases and decreases in input temperature to have the same effect on ER mod_surf . Finally, for the orthomosaic from July 2018 in Fig. 8a, the temperature error decreased ER mod_surf . Since ER mod_surf peaked at ~40 • C in 2018 (Fig. 5a), which was also the mean orthomosaic surface temperature in Fig. 8a, both increases and decreases in input temperatures of individual pixels caused declines in ER mod_surf . The effect of the temperature measurement uncertainty on upscaled ER mod_surf is thus closely linked to the mean surface temperature and the shape of the ER model used for the upscaling. We also quantified the impact of the spatial distribution of the two vegetation communities when upscaling ER by calculating the mean ER mod_surf for each orthomosaic (Figs. 8 and 9) based on 1) the land cover classification and separate ER model fits for hummocks and hollows (i.e. ER mod_surf ); and 2) a single ER model fit (for both vegetation communities combined, 'ER mod_surf_all ', Fig. 5g-h) for the whole orthomosaic (Table S7). All the orthomosaics had a much larger proportion of hummocks (up to 85% coverage) than hollows (up to 21% coverage; Table S7). The comparison showed that neglecting the ER differences between the two vegetation communities caused mean ER mod_surf to change by up to ±18% (Table S7). ER mod_surf_all both overestimated (2018) or underestimated (2019) the true ER mod_surf, depending on whether hollows or hummocks had higher ER.

The choice of temperature metric affects ER mod temperature sensitivity
The most significant impact of using different temperature metrics to model ER was on the temperature-sensitivity of ER mod . Several studies have shown that the temperature sensitivity of soil respiration or ER increases with measurement depth when using soil temperature, or when using soil as opposed to air temperature (Graf et al., 2008;Lasslop et al., 2010;Subke and Bahn, 2010), including in peatlands (Lafleur et al., 2005). This effect is caused by the decreasing range of soil temperatures with increasing soil depth. The differences in temperature sensitivity can affect the propagation of uncertainties in ER mod . For example, ER mod_soil had a very high temperature sensitivity, such that even a small shift in input temperature would cause a large difference in the estimated flux. Such shifts may be important for sites that do not have uniform soil temperatures and only have a few measurement points across the site. ER mod_surf , on the other hand, had the lowest temperature sensitivity, which is expected given the greater temperature amplitude experienced at the peat or plant surface compared to in the air or soil (see Section 3.3). However, surface temperature measurements currently have higher uncertainties than air or soil temperature measurements (see Sections 4.2 and 4.3 ), and further technological improvements would hence highly benefit this model approach. The choice of temperature metric also impacted estimates of mean growing season ER mod (for the collars), which was up to 16% higher when air (compared to surface) temperature was used to model ER.

Surface temperature yields similar ER mod accuracy to conventional temperature metrics
We found similar ER mod accuracy when using surface temperature compared to soil or air temperature to drive the model. Furthermore, when we accounted for the temperature measurement uncertainties in the thermal camera data, mean growing season ER mod_surf changed by <5%. Thermal camera data thus provide an alternative to conventional temperature metrics for ER modelling when high spatial resolution temperature information is desired (see Section 4.3). These findings are reinforced by Pau et al. (2018) who found that using thermal imagery-derived canopy surface temperatures improved GPP model accuracy compared to using air temperature. Our results also support the use of airborne and satellite land surface temperature data for upscaling ER mod_surf at larger scales (though see Section 4.7 regarding the impact of drought on surface temperature).

ER mod_surf can be upscaled in unprecedented detail and account for spatial variations in ER
Using thermal data to model ER offers new opportunities to upscale and map how ER mod_surf varies across an ecosystem, because thermal images can be collected with a UAV at very high resolution (<7 cm). The importance of accounting for the spatial variability of ER across an Fig. 9. Upscaled ER mod_surf using UAV imagery, showing diurnal changes in temperature and ER mod_surf . Thermal orthomosaics of surface temperature (a-c), with the date and starting time of each flight. ER mod_surf maps (d-f) using temperatures from (a-c) as input to the hummock and hollow GPM model fits for 2019. ecosystem was highlighted by our findings that hollows and hummocks have significantly different ER and also responded differently to the hot drought. Our analysis of the UAV-derived ER maps showed that mean ER mod_surf was over-or underestimated by up to 18% if this ER difference was not accounted for. Using a single parameter set to model or upscale ER would mask the differing reactions of vegetation communities to changes in abiotic variables or extreme events like droughts and hot spells. The significance of such biases in other peatlands will depend on the relative areal coverage of each vegetation type and the differences in their ER. Having high resolution vegetation classifications derived from e.g. UAV imagery, as well as separate model fits for each vegetation type, is therefore crucial to enable accurate upscaling.
High resolution ER maps provide a novel resource to examine the key drivers of the spatial variability of ER within the footprint of eddy covariance flux tower measurements, and can do so at a much higher spatial resolution than e.g. chamber measurements. In addition, ER maps could be used to assess whether flux tower measurements are representative of the surrounding ecosystem and to provide estimates of ER over larger areas (e.g. a whole peatland) than can be covered using flux tower measurements. Such information can help minimize or explain biases introduced by a mismatch of coarse spatial resolution satellite data and the tower flux footprint (Gelybó et al., 2013).

Accounting for the impact of drought on ER
We observed that the ER of both vegetation communities decreased during the 2018 hot drought. Previous research generally suggests that peatland ER increases during warmer and dryer periods due to the positive relationship between ER and temperature, as well as the increasing availability of oxygen for soil respiration as the water table drops (Dorrepaal et al., 2009;Fenner and Freeman, 2011). Nevertheless, there is also evidence that extreme water-table drawdown or an extended dry period can lead to declines in peatland ER (Laiho, 2006;Lund et al., 2012;Strack et al., 2009). The 2018 hot drought was an extreme event, both in terms of its longevity and the combination of high air temperatures and low precipitation at the peak of the growing season, leading to an extended period of water table drawdown (Peters et al., 2020;Rinne et al., 2020). Easily decomposable labile C exposed by the lower water table may have been used up in the early stages of the growing season, leaving only the more recalcitrant fraction of organic matter and thus slowing decomposition rates (Laiho, 2006;Strack et al., 2009). Furthermore, both heterotrophic and autotrophic respiration decline under water stress (Atkin and Macherel, 2009;Moyano et al., 2013). Hummocks would have been more affected by reductions in autotrophic respiration due the greater plant biomass than in hollows. A decline in plant productivity due to the extremely dry conditions (Goodrich et al., 2015;Otieno et al., 2012) would also have reduced the vascular plant 'priming' effect on soil microbial respiration (Gavazov et al., 2018;Walker et al., 2016), thus further reducing hummock ER. Finally, ER rates may have acclimated to the consistently high temperatures, which would lead to a decrease in the temperature sensitivity of ER and thus to reductions in ER (Atkin and Tjoelker, 2003).
We note the importance of the function shape used to model ER for predicting whether there will be rapid increases in ER during extreme drought and heatwave events. We found that, across all three temperature metrics, ER stopped increasing at high temperatures during the hot drought (with the only exception of ER mod_soil for hollows, Fig. 5). We chose to use the GPM model because it captured the observed plateauing of ER at high temperatures during the hot drought, whereas other, more commonly-used models (e.g. Lloyd and Taylor, 1994) assume continuous increases in ER with temperature. Decisions around which model to use are particularly important given that the temperature-ER relationship is commonly used to partition and gapfill flux tower data. Using models that assume continuous increases in ER at high temperatures may lead to overestimates of carbon loss during severe droughts when water-stress can limit ER.

Surface temperature provides insights into the impact of drought
Canopy surface temperature is linked to stomatal conductance and evapotranspiration because these processes affect heat exchange at the leaf surface. As a result, surface temperature has been used to quantify plant drought stress, particularly in precision agriculture (Berni et al., 2009;Maes and Steppe, 2012;Still et al., 2019). We observed significant differences in surface and soil temperature between the hummocks and hollows during the 2018 hot drought. During these periods, soil temperature was warmest in hummocks, but surface temperature was >3 • C warmer in hollows than hummocks. These opposite changes in soil and surface temperature could be evidence of a negative feedback mechanism to minimize peat water loss (Kettridge and Waddington, 2014). Under this mechanism, the surface resistance of moss increases during drought, thus minimizing water loss by evaporation and desiccation of the underlying peat. The increased surface resistance is associated with a rise in surface temperature due to the decreased latent heat flux. Thus hollows may have had very high surface temperatures but maintained a more stable water table level (and lower soil temperature) whereas vegetation on hummocks continued transpiring (lower surface temperature) which led to rapid decreases in water table level (Fig. 3c) and higher soil temperatures.
This negative feedback mechanism may also explain differences in hollow ER mod_surf and ER mod_soil during the hot drought. The decoupling of surface and soil temperatures caused surface temperatures to continue increasing to almost 50 • C, whilst increases in soil temperatures were dampened, reaching a maximum of 28 • C. Thus hollow ER mod_surf plateaued at high surface temperatures, whereas hollow ER mod_soil showed sustained increases in ER with soil temperature. Using thermal images thus enables a more comprehensive analysis of how an ecosystem responds to drought and the effects of drought on water and energy fluxes.

Limitations and future work
We presented the results of a pilot project to assess the feasibility of using thermal cameras and UAVs to model and upscale the ER of peatlands. To further develop our modelling and upscaling approach, several limitations must be overcome which we discuss below and use to highlight opportunities for further work.
Thermal data, particularly from UAVs, should be used cautiously, since deriving accurate temperature data from these cameras is not trivial (Kelly et al., 2019). Our uncertainty analysis showed that the inclusion of a ±5 • C random error changed mean ER mod_surf of the orthomosaics by ±7%, depending on the mean surface temperature and shape of the ER model used for the upscaling. The uncertainties in the thermal data could be minimized by generating custom calibration equations for the camera (Ribeiro-Gomes et al., 2017) or by correcting for temperature drift (Mesas-Carrascosa et al., 2018). These additional steps may be avoided by further technological improvements of thermal cameras and by using camera models with higher accuracy specifications.
The thermal data processing is also affected by the choice of emissivity value: a 1% change in emissivity can lead to a ~1 • C change in measured temperature (Aubrecht et al., 2016). Although emissivity varies depending on the density, structure, water content and species of the vegetation, using a fixed value, as applied in the present study, is often necessary due to the challenges of measuring the emissivity spectra of vegetation, particularly under field conditions (Olioso et al., 2007;Ribeiro da Luz and Crowley, 2007). Measurements of emissivity spectra for peatland vegetation are rare, and experimental data on the effect of water stress have shown both increases and decreases in vegetation emissivity of ~0.02 (Buitrago et al., 2016;Olioso et al., 2007). We estimate that the use of a fixed emissivity value (0.97) may have caused an error of roughly ±2 • C in our surface temperature measurements (assuming an emissivity value of 0.99 for flooded areas and a decline in emissivity of 0.02 due to water stress). Nevertheless, it is likely that the large uncertainty range associated with the thermal cameras used in this study (±5 • C), as well as limitations in the land cover classification and ER model (discussed below) would have overshadowed errors introduced by the choice of emissivity value.
As discussed in Section 4.3, accounting for the spatial distribution of hollows and hummocks is key to accurately upscaling of ER, and thus the land cover classification used has an important impact on the upscaling results. Becker et al. (2008) suggest a minimum spatial resolution of 60 cm to accurately upscale CO 2 fluxes from peatland microforms, which we surpassed with our land cover classification data (3.5 cm resolution). The accuracy of our classification was also similar to their work and that of Lehmann et al. (2016), who used an object-based classification of peatland vegetation. However, the confusion matrix (Table S2) suggests that the hollow class was incorrectly classified as hummock for 32% of the validation points used (but the opposite error occurred in only 1% of cases), indicating that the proportion of hollows in the peatland may be underestimated. Therefore, the ER mod fluxes in Figs. 8 and 9 may have been under-and overestimated in 2018 and 2019, respectively. This issue may be resolved in future work by performing a land cover classification with a larger number of vegetation or microtopographical classes, further improving the accuracy of the modelled and upscaled greenhouse gas fluxes (Moore et al., 2019).
We note that the accuracy of the ER mod maps presented in Figs. 8 and 9 is also dependent on the robustness of the ER measurements and accuracy of the ER models used for the upscaling. Studies based on chamber measurements can face challenges due to the spatial and temporal autocorrelation of the measurements because collar placement and measurement timing are often constrained by logistics (e.g. site accessibility, power availability, location of boardwalks) and this study is no exception. The collars used for the manual ER measurements were clustered in the field of view of the thermal camera to allow a direct link between the surface temperature and ER measurements. Our upscaling method assumes that these measurements are representative of the whole site. The temporal autocorrelation of chamber greenhouse gas measurements is rarely considered. Nevertheless, future modelling approaches could use mixed effects models to account for these issues (Kravchenko 2015).
In terms of model accuracy, our results for ER mod_surf (mean R 2 = 0.54, mean NMRSE = 47%) are in line with previous work using chamber data in peatland sites and simple empirical ER models based on temperature: R 2 between 0.2-0.9 and NRMSE between 15-50% (Acosta et al., 2017;D'Angelo et al., 2016;Juszczak et al., 2013;Mäkiranta et al., 2009). We modelled ER using only temperature data, which is a common approach in the greenhouse gas flux and remote sensing communities (Ai et al., 2018;Reichstein et al., 2005;Reichstein and Beer, 2008). However, the goodness-of-fit statistics for our model fits indicate that more variables could be included to better capture the variability in the observed ER. Drought clearly impacts ER and we found strong relationships between ER and water table depth (particularly in hollows). UAV data has recently been used to estimate soil moisture content through the Temperature Vegetation Dryness Index method (Wigmore et al., 2019) as well as water table depth based on point-cloud and land cover classification data (Rahman et al., 2017). ER models including soil moisture are available (e.g. Carlyle, 1988;Tang et al., 2005) and these would eliminate the need to parameterize the ER model separately for drought and non-drought years as we have done here. Furthermore, ER is closely linked to vegetation productivity which can be related to spectral indices representing vegetation greenness (Jägermeyr et al., 2014). Combining water table depth, vegetation indices and surface temperature data from ground-based sensors and from UAV flights is a logical next step in improving our ER modelling and upscaling approach.

Conclusions
In this study, we present a novel approach for upscaling ecosystem respiration (ER) at a hemi-boreal peatland in Sweden. The approach is based on surface temperature measurements from a tower-based thermal camera in tandem with a UAV-mounted thermal camera, and offers maps of ER at unprecedented high spatial resolution (<7 cm). The suitability of surface temperature as a driver for modelling peatland ER was successfully validated by a comparison with soil and air temperature-driven ER models at the chamber scale.
The resulting maps reveal how ER can vary significantly across an ecosystem at the spatial scale of the dominant vegetation communities. Our findings hence emphasize the importance of modeling ER separately for the two vegetation communities present on the peatland. The need to account for vegetation heterogeneity at the site was further highlighted by the different response of the two communities to the extremely hot and dry summer in 2018. During this extended hot drought, ER declined in both communities but more in hummocks than in hollows. These results suggest that peatland carbon loss during severe drought may, in some cases, be lower than expected, and that it will be highly dependent on vegetation composition.
The present study demonstrates the potential of surface temperature data from field-based thermal cameras as a new tool for modelling and upscaling the ER of peatlands. Future work may include productivity and soil moisture information in remote sensing ER models, to improve the models' ability to capture how the spatial heterogeneity of ER affects the aggregated response of a peatland to extreme weather events and climate change. The robustness of this method will be further enhanced through future technological improvements to the accuracy of thermal camera measurements. Ultimately, our upscaling approach lays the foundation for a more detailed investigation of how ER varies across an ecosystem, to assess the representativeness of flux-tower derived ER by comparing vegetation-dependent ER within its flux footprint with the ecosystem estimate, and to interpret biases that occur when using coarse resolution satellite data to upscale tower-based measurements.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.