High-resolution analysis of observed thermal growing season variability over northern Europe

Strong historical and predicted future warming over high-latitudes prompt significant effects on agricultural and forest ecosystems. Thus, there is an urgent need for spatially-detailed information of current thermal growing season (GS) conditions and their past changes. Here, we deployed a large network of weather stations, high-resolution geospatial environmental data and semi-parametric regression to model the spatial variation in multiple GS variables (i.e. beginning, end, length, degree day sum [GDDS, base temperature + 5 °C]) and their intra-annual variability and temporal trends in respect to geographical location, topography, water and forest cover, and urban land use variables over northern Europe. Our analyses revealed substantial spatial variability in average GS conditions (1990–2019) and consistent temporal trends (1950–2019). We showed that there have been significant changes in thermal GS towards earlier beginnings (on average 15 days over the study period), increased length (23 days) and GDDS (287 °C days). By using a spatial interpolation of weather station data to a regular grid we predicted current GS conditions at high resolution (100 m × 100 m) and with high accuracy (correlation ≥ 0.92 between observed and predicted mean GS values), whereas spatial variation in temporal trends and interannual variability were more demanding to predict. The spatial variation in GS variables was mostly driven by latitudinal and elevational gradients, albeit they were constrained by local scale variables. The proximity of sea and lakes, and high forest cover suppressed temporal trends and inter-annual variability potentially indicating local climate buffering. The produced high-resolution datasets showcased the diversity in thermal GS conditions and impacts of climate change over northern Europe. They are valuable in various forest management and ecosystem applications, and in adaptation to climate change.


Introduction
Climate change has led to major alterations in global air temperature and precipitation patterns with asymmetric manifestation between seasons, and these changes are projected to further continue over the upcoming decades (IPCC 2013; Bintanja and Andry 2017). As a consequence, changes in thermal growing season (GS), that is, the period of suitable conditions for plant growth, have been widely observed and projected over northern hemisphere (Linderholm et al. 2008;Jeong et al. 2011;Ruosteenoja et al. 2016Ruosteenoja et al. , 2019. In general, such changes have included earlier beginning and later end of GS, with consequent increases in the length of GS and accumulated growing degree day sum (°C days, GDDS). These changes in thermal GS will profoundly impact on ecosystem functions, biodiversity patterns, agricultural and forest management practices, bioeconomical activities as well as exposure to natural hazards (Venäläinen et al. 2001;Piao et al. 2007;Linderholm 2006;Liu et al. 2018;Jactel et al. 2019). Consequently, there is an urgent need for up-todate and spatially detailed information of current variability and recent changes of thermal growing season conditions. Northern Europe experiences large spatial and temporal variability in weather and climate due to its position in the end of North Atlantic storm track (i.e. the Polar Front), Scandes mountain, proximity to Atlantic ocean in the west and large Eurasian continent in the east (Tikkanen 2005;Wernli and Schwierz 2006). Moreover, its distinct snow climate and specialized ecosystems makes it especially vulnerable to climate change (Aalto et al. 2017a;Niittynen et al. 2018;Niskanen et al. 2019). Boreal environments are in general humid and moist, particularly in the spring after snow melt. Thus, water availability is not limiting the onset of the thermal growing season. The northern Europe (i.e. Finland, Sweden, Norway and Denmark) is a unique region in that it is covered by a dense and reliable network of weather observation stations operated already in the 1950s, that is, before the main satellite era. Such a well-established observational network is a prerequisite for accurate tracking of climate change and conducting detailed spatial analyses of climate variability (Tveito et al. 2001;Aalto et al. 2017b). Weather stations are often placed on open landscape to represent broad environmental conditions (De Frenne and Verheyen 2016). Recent studies show that, with sufficiently large station network, such data may be used to model local air temperature, and further assess thermal growing season variability in respect to local conditions such as topography and water cover variables over large spatial domains (Meineri and Hylander 2016; Aalto et al. 2017b).
Current knowledge of GS variability over northern Europe is founded on partly outdated analyses, short analysis periods, limited geographical coverage and/or coarse spatial resolution. For example, Tveito et al. (2001) presented a collection of climate maps over Nordic areas including also beginning (hereafter GS beg ), end (GS end ) and length (GS len ) of GS over the period of . Surprisingly, no updated GS maps or analyses of past trends covering the northern Europe have been published since [but see Linderholm et al. (2008) for the Greater Baltic Area]. This knowledge-gap limits our understanding of the current pace of climate change and its potential ecosystem impacts over northern Europe. In more thorough but spatially limited analysis, Irannezhad and Kløve (2015) used gridded daily weather data to quantify changes in GS beg , GS end and GS len over Finland 1961Finland -2011. The authors related the spatial variability in temporal GS trends to various atmospheric teleconnection patterns. They found significant correlations with multiple oscillation indices indicating that interannual GS variability is responsive to large-scale atmospheric circulation dynamics. However, the use of rather coarsescale climate data (spatial resolution of 10 km × 10 km) masks out much of the landscape level GS variability important for understanding changes that local ecosystems are facing (Aalto et al. 2017b). In addition to thermally defined GS, vegetation phenology-based GS estimates have been produced over northern Europe (Karlsen et al. 2007(Karlsen et al. , 2008(Karlsen et al. , 2009Høgda et al. 2013). The use of remote sensing techniques allows for improving the spatial details of the analyses (e.g. 250 m for Moderate Resolution Imaging Spectroradiometer [MODIS]), but so far such investigations have been limited to relatively short periods of past 20 years hindering understanding of long-term GS dynamics. Therefore, there is a clear gap in understanding current and past GS variability in Northern Europe.
In general, GS variability follows air temperature patterns that in turn can show substantial local scale variability depending on e.g. topographical settings and proximity to water bodies (Ashcroft and Gollan 2012;Meineri and Hylander 2016;Aalto et al. 2017b). In addition to strong effects of geographical location (i.e. latitudinal and longitudinal position) and elevation gradient on air temperatures via lapse rate (Rolland 2003), local scale variation in incoming solar radiation intercepted by topographical and forest canopy shading together with effects due to e.g. latent heat fluxes can create substantial spatial heterogeneity in GS variables (Fridley 2009;Dobrowski 2011;Meineri and Hylander 2016;Greiser et al. 2018;De Frenne et al. 2021). Lakes are prominent in the northern Europe due to multiple glaciations and they are known to influence local climate by suppressing temperature variations due to their high specific heat capacity that leads to a cooling effect during early GS (i.e. delaying beginning) and warming effect during late GS (delaying end, Lookingbill and Urban 2003). Such local effects may lead to temperature and GS trends that by their magnitude and spatio-temporal dynamics can differ from regional patterns, with potentially important consequences on local ecosystem dynamics (Pepin and Seidel 2005;Daly et al. 2010;Dobrowski 2011;Lenoir et al. 2016). Moreover, studies suggest that the Urban Heat Island effect (UHI; elevated urban temperatures compared to surrounding countryside; Oke 1973) can exert a control on thermal growing seasons and phenology in urban areas, where urban surfaces and structures alter energy balance, heat fluxes and air-mixing (Oke 1995;Jochner et al. 2012;Zipper et al. 2016). However, contemporary gridded GS data does not capture such local scale heterogeneity impeding our ability to understand GS dynamics and their ecosystem impacts in northern Europe.
Here, we document and analyze current (1990-2019) thermal growing season conditions, their interannual variability and past changes (1950-2019) over northern Europe using a large network of weather stations. We provide the most up-to-date and spatially detailed picture of the current and past variability in multiple GS variables (i.e. GS beg , GS end , GS len and GDDS) by statistically modeling their spatial variation in respect to geographical location, topographical, water and forest cover, and urban land use variables. Finally, we produce geospatial GS data layers covering northern Europe at very high spatial resolution (100 m × 100 m) to be used in various ecosystem applications.

Study domain
We analysed GS variability in northern Europe, within the domain of ca. 55-71° N, 5-35° E (Fig. 1). The climate of this region is dominated by the effects of proximate Atlantic Ocean, annual and decadal variability of North Atlantic Oscillation, moving low pressure systems at the Polar front, elevation due to Scandes Mountain range and large Eurasian continent (Hurrell 1995;Tikkanen 2005). Based on the daily gridded climate dataset (Nordic Gridded Climate Dataset [Tveito et al. 2005, updated], spatial resolution of 1 km × 1 km, https:// surfo bs. clima te. coper nicus. eu/ dataa ccess/ access_ ngcd. php), the average annual air temperature conditions (1990-2019) vary from − 7.8 to 9.6 °C. The study area has extensive latitudinal and elevational gradients leading large variations in incoming solar radiation. Due to recent glaciations, lakes are prominent in the study area.

Weather station data
Weather station data for Finland, Sweden, Norway, Denmark and westernmost parts of Russia were compiled from the European Climate Assessment & Dataset (ECA&D, accessed 19.6.2019, total of 1395 stations) data base (Klok and Klein Tank 2009) and represent daily mean air temperature (Tday) data 1950-2019 ("non-blended" data series). The data have passed national operative quality control schemes. Russian stations were included to the dataset to reduce border effect in subsequent spatial modeling as well as to increase the environmental representativeness of the station data.
The spatial distribution of the observations is dense enough for studying daily aggregated values, especially when those daily values are used to compute annual quantities (average distance < 200 km; Fig. 1). Such station density is adequate in characterizing weather phenomena of the region having the smallest spatial scales between 20 and 200 km (Thunis and Bornstein 1996). Despite this, there are regions, such as the mountains in South Norway, where the network is sparser and where greater uncertainty in the reconstruction of the gridded fields can be expected.

Response variables: growing season variables
We analyzed four variables that captures the key properties of thermal growing seasons: beginning of the GS (GS beg , day of the year [DOY]), end of the GS (GS end , DOY), length of the GS (GS len , days) and growing degree day sums (GDDS, °C days). We define the thermal growing season as the period when daily mean temperature (Tday) is permanently at or above + 5 °C (Eqs. 1-2). The beginning and end of the growing season was determined using the so-called integral method (see Ruosteenoja et al. 2016), which identifies the date after the absolute minimum of the sum(Tday-threshold) has been reached (GS beg ) and analogously GS end when the absolute maximum of the sum(Tday-threshold) has been reached, but not earlier than 1st of June.
A minimum of 95% of the annual Tday data must exist in order to merit further calculations of GS variables.
To characterize the current spatial variability in GS conditions over the study area, we averaged the four annual GS variables over the period of 1990-2019 (hereafter denoted as µGS beg , µGS end , µGS len and µGDDS). A minimum of 10 years of data was required to calculate the average values. This liberal limit was chosen in order to increase stations' coverage over the environmental gradients especially towards high elevations (Fig. S1), where the observational network is sparser. Further, this preselection yielded 482 stations for consequent analyses with an average of 26 years of data (Fig. 1b).
We estimated the temporal trends in GS conditions (βGS beg , βGS end , βGS len and βGDDS) over the period of 1950-2019 using Sen's slope estimator (Sen 1968) (function sens.slope in R package trend) (Pohlert 2020) for stations having at least 45 years of data. This preselection resulted in trend estimates for 349 stations to be used in consequent spatial analyses, which on average have 61 years of data ( Fig. 1c), and on average 28 years of data covering the three most recent decades (i.e. 1991-2019). Using the same subset of 349 stations, interannual variability of the GS variables (σGS beg , σGS end , σGS len and σGDDS) were quantified in terms of standard deviation over the period of 1950-2019.

Predictors: geospatial data
We used several geospatial data to derive predictor variables (hereafter referred as predictors) for explaining the spatial variation in the response variables. The geospatial data were managed using ArcMap (version 10.7.1.) Spatial Analyst-functions.
Topography exerts a strong control on air temperature and thus GS variability via elevation (temperature lapse rate) and solar radiation interception (surface energy input). For Fig. 1 Weather station network. a Weather stations recording daily air temperatures that were used to calculate temporal trends (β) and interannual variability (σ, standard deviation; 1950-2019, red squares, n = 349) and mean conditions over 1990-2019 (µ, black squares, n = 482 including also the previously mentioned 349 stations) in thermal growing season variables. The background map depicts the elevation above the sea level (m a.s.l.) and grey coloring indicates areas outside the analysis domain. b-c Histograms showing the amount of years available for calculating µ (b), and β and σ (c), respectively. d-g Examples of time series of thermal growing degree day sum (GDDS), beginning of the growing season (GS beg ), end of the growing season (GS end ) and length of the growing season (GS len ), respectively, from selected weather stations (locations indicated in a). DOY day of the year. Statistical significance of the temporal trends (β, estimated using Sen's slope) is indicated as *** = p ≤ 0.001, * = p ≤ 0.05, n.s. = not significant (p > 0.05). The black lines in d-g depict the average of 1990-2019 and the blue lines depict the temporal trends fitted using spline smoothing obtaining elevation information we used the values reported in the station metadata. Other topography-related predictors were derived from the Merit DEM (Multi-Error-Removed Improved-Terrain Digital Elevation Model) that covers land areas between 90° N and 60° S at the spatial resolution of 90 m × 90 m (Yamazaki et al. 2017). We calculated several potential predictors from the DEM that have previously been used to explain spatial variation in local air temperatures at the study area (Meineri and Hylander 2016;Aalto et al. 2017b). To capture the effect of solar radiation on air temperatures, we calculated the northness index, which is a cosine-transformation of an aspect layer (in degrees). The resulting layer has values ranging from − 1 (south-facing slope) to 1 (north-facing slope). We also calculated potential incoming solar radiation (PISR, MJ cm −2 a −1 ) using the equations in (McCune and Keon 2002). However, this potential predictor was excluded from the analyses due to very high positive correlation with latitude that confounds the multivariate models (Meineri and Hylander 2016). We also considered the effect of stations' local topographical position as a potential predictor for capturing the effect of local topography on air-mixing (Ashcroft and Gollan 2013). Following the methodology by Aalto et al. (2017b) we calculated relative elevation as the difference of a cell to the minimum elevation at the radius of 500 m and 2000 m. Here, we assumed that locations that have relative elevation close to zero are prone to more stable atmospheric conditions than locations above the local minimum (Pepin et al. 2009). However, after preliminary testing this potential predictor (both with 500 m and 2000 m radius) was omitted for further analyses due to highly positively skewed distribution (results not shown).
Air temperatures in the study area are strongly affected by the proximity of the North Atlantic Ocean, The Baltic Sea and the Barents Sea (Tikkanen 2005). Following the methodology in Meineri and Hylander (2016) we calculated distance to sea (in meters, using ArcGIS Spatial Analystfunction Cost Distance), and to simplify our approach, we assumed that all surrounding sea areas have similar effect on GS variables. However, due to varying ocean dynamics e.g. in freeze-up and circulation this assumption may not hold and can introduce additional uncertainty to our analyses. We used global water cover data (World Water Bodies; UCLA Institute for Digital Research and Education 2019) to delineate lakes and shorelines. Following the methodology in Aalto et al. (2016), Meineri and Hylander (2016) and Aalto et al. (2017b) we calculated distance to lakes (lakes ≥ 10 km 2 considered) to capture local effects of lakes on daily air temperatures and further GS variability (Lookingbill and Urban 2003). To account for the effects of forest cover on air temperatures via mediating net radiation, air flow and latent heat fluxes, we used the global data product by Hansen et al. (2013) which represents a remotely-sensed estimate of the percentage forest cover (0-100%) in 2000 at the spatial resolution of 30 m × 30 m. Finally, we used the data by Gao and O'Neill (2020) to depict the fraction of urban land to all land areas (0-1) in 2000 (https:// datav erse. harva rd. edu/ file. xhtml? persi stent Id= doi: 10. 7910/ DVN/ ZHMI1L/ TEVEB 0& versi on=1.0; accessed 9.6.2021). The spatial resolution of the urban fraction data is 0.125° × 0.125° (ca. 12.5 km at the equator). Despite this relatively coarse spatial resolution, we find the data suitable for our modeling in quantifying the general effect of urban land use on air temperatures and consequently on thermal GS. For the subsequent modeling, the forest cover and urban fraction variables were log(x + 1) and log(x + 0.01) transformed to balance their highly positivelyskewed distributions (Fig. S1), respectively. In addition to predictors described above, we used geographical location (latitude and longitude) to describe the regional-scale variations in GS variables caused by e.g. maritime-continental gradient, seasonal variation in cloud cover and broad scale advection affecting GS temperature patterns.
Water-related predictors were rasterized, and all geospatial data were reprojected to same coordinate reference system ETRS89-LAEA (epsg: 3035) and further resampled to a common spatial resolution of 100 m × 100 m using bilinear interpolation (ArcGIS resample-function).

Statistical modeling
All statistical analyses were conducted using the R software version 3.6.1. (R Development Core Team 2011). The response variables (i.e. µGS, βGS and σGS) were related to the predictors using generalized additive modeling (GAM; Hastie and Tibshirani 1990) as implemented in R-package mgcv (Wood 2011). The GAM are semi-parametric extensions of generalized linear models that use smoothing functions to fit non-linear response functions to the data: where g( ) is the link function that connects the estimated mean to the distribution of the response variable (here we assume Gaussian errors), 0 is the intercept, s i is the smoothing function to be estimated and x i is a predictor.
Following Meineri and Hylander (2016) we used a twostaged modeling procedure: first, the response variables were modeled as a function of topographical, water and land cover predictors (Model1, Eq. 4), and then the residuals of model1 (e) were modeled as a function of geographical location (Model2, Eq. 5) using anisotropic tensor interaction term (te-function) that allow for an asymmetric effect of geographical location (Wood 2011): Here k refers to the maximum smoothing function which is further optimized by the model fitting algorithm (Wood 2011). To reduce potential overfitting the k-argument was set to three. The two-staged modeling procedure was done so that topographical, water and land cover predictors can contribute in explaining the variation in the response variables, which would otherwise be masked out by the strong effect of geographical location on the GS variables. We did not use a model selection criterion i.e. all considered predictors were kept in the models regardless of their statistical significance or other optimization criteria. This was done because (i) the most of the chosen predictors have previously been used in spatial modeling of air temperature conditions in the study region (Aalto et al. 2016(Aalto et al. , 2017bMeineri and Hylander 2016), and (ii) the estimated effects of local-scale predictors, especially northness and forest cover, are expected to be modest due to weather stations' placement on open and relatively flat landscapes. As a supplementary analysis of spatial variation in temporal trends (βGS), we included the number of years as an additional predictor to the model. This was done in order to test whether this feature of the data (i.e. varying amount of years for trend calculations) influences the observed spatial patterns of temporal GS trends.
Effect sizes (EF) for each predictor, after holding other predictors constant at their mean values, were calculated by subtracting the predicted maximum GS values from the predicted minimum GS values (Aalto et al. 2017b). Thus, EF represents the magnitude of variation in a response variable caused by a predictor. Finally, the models were used to predict the GS variability over the study area at the spatial resolution of 100 m × 100 m. For each response variable (µGS, βGS and σGS) this was done by summing the predictions of the Model1 and the Model2 into final predictions. To secure consistency of the final data layers, µGS len was produced by subtracting the predicted µGS end from the predicted µGS beg .
To assess the accuracy of the spatial predictions we used leave-one-station-out cross-validation scheme, where each station in turn was left aside from the observation data, the model was fitted using n − 1 observations and consequently predicted over the station withheld from model fitting. This (4) GSvariable = f (elevation, northness, distance to sea, distance to lake, forest cover, urban fraction) procedure was repeated as many times as there were stations in the data. The prediction accuracy was assessed in terms of Pearson's correlation coefficient (r), mean error (bias) and root mean squared error (rmse) between observed and predicted values (Aalto et al. 2016).
Spatial autocorrelation (SAC) is a common property of any spatial dataset and means that observations are related to one another by the geographical distance (Legendre et al. 2002). SAC in the model residuals violates the independence assumption commonly required by statistical models and may lead to inflated hypothesis testing and biased model estimates (Beale et al. 2010). To investigate whether GS observations and model residuals were spatially autocorrelated, we calculated correlograms (as implemented in R-package pgirmess version 1.6.9. (Giraudoux 2018)) which describes the spatial dependency between the observations as a function of distance between the point pairs in terms of Moran's I.

Prediction accuracy of the spatial models
The tested predictors were only moderately inter-correlated (pairwise Pearson's correlation coefficient ≤ 0.55) and thus our models were not confounded by multicollinearity ( Fig.  S3-S4). The statistical models were able to predict µGS with high accuracy (Fig. 3a-d) with correlation (r) between observed and predicted values ≥ 0.92. Spatial variation in βGS was more challenging to predict as for example the model for βGS end resulted in poor agreement with observations (r = 0.33; Fig. 3f). The other βGS predictions showed reasonable agreement with correlations ranging from 0.56 (βGS len ; Fig. 3g) to 0.63 (βGS beg and βGDDS; Fig. 3e, h, respectively). Prediction accuracy for interannual variability σGS was reasonable with correlations ranging from 0.46 (σGS end ) to 0.76 (σGS beg ) (Fig. 3i-l).

Spatial GS patterns and their drivers
In general, spatial variation in µGS was strongly constrained by geographical location with clear latitudinal gradients and elevation (Figs. 4, 5). However, the predictions suggest substantial local variability in µGS (Fig. 6). For µGS beg we found that northness has a positive effect (EF = 1.50 days) after holding other predictors at their mean values. On average µGS beg increased along an increasing distance to sea (EF = 8.69 days) whereas for µGS end and µGS len we found an opposite relationship. Both µGS end and µGS len were found to be positively related to the distance to lakes (EF = 27.51 days and 35.74 days, respectively) and negatively to forest cover (EF = − 2.86 days and − 3.96 days, respectively). Urban fraction was found to be negatively related to µGS beg (EF = − 17.25 days) and positively related to µGS end (EF = 12.71 days) and µGDDS (EF = 366.61 °C days).
Our analyses indicated that on average βGS beg has decreased (i.e. earlier beginning) the most in southern parts of the study area (over five days per decade) and at low elevations (Figs. 7, 8, Fig. S2). Moreover, the modeling suggests that βGS beg is positively related to distance to sea (EF = 1.45, when trends are expressed as per decade; Fig. 8c) while increasing distance to lakes tends to lead to larger  Fig. 8d) . While there are considerable uncertainties in the modeled spatial patterns of βGS end (Fig. 3f), our modeling identified positive links (i.e. increasing trend) with elevation (EF = 0.54), distance to sea (EF = 0.66) and distance to lakes having the largest effect (EF = 1.34). βGS end was found to be negatively related to forest cover and urban fraction. The spatial patterns of βGS len are mostly constrained by geographical location (EF = 3.17) and a positive effect of distance to lakes (EF = 2.65). βGDDS showed clear negative trends along latitudinal and elevational gradients having the largest effect. Further, our results suggest that βGDDS is negatively related to the distance of the sea (EF = − 7.89), forest cover (EF = − 5.65) and positively related to the distance of lakes (EF = 22.15) and urban fraction (EF = 8.56). In supplementary analysis we found that the number of years of data were significantly related to all four βGS variables (Fig. S5) indicating smaller temporal trends over stations with more data. However, the estimated effects (response shapes and effect sizes) are in good agreement with the model results neglecting this predictor (Fig. 8) and thus we consider our initial modeling being valid.
We found that in general the spatial variation of σGS was related to elevation gradient (negative effects for σGS end , σGS len and σGDDS) and distance to sea (negative effect) as well as distance to lakes (positive) (Figs. S6-S7). All GS variables and the residuals of Model1 (Eq. 3) showed significant spatial autocorrelation patterns (Fig. S8). It is also worth remarking Fig. 5 (panels a, g, m, and s) shows that the uncertainty of the GS variables increases at higher elevations, most likely because of the sparser observational network available there, as pointed out in the previous sections.

Comparison of the estimated growing season trends with previous studies
Our results highlight the remarkable spectrum of growing season conditions over northern Europe and consistent temporal trends (e.g. GDDS increase ranging 84-602 °C days over the study period of 1950-2019) that may be mediated by several environmental gradients. In agreement with previous research from the area, we found that the increase in GS len can mainly be attributable to earlier beginning of GS (Linderholm et al. 2008;Irannezhad and Kløve 2015), as we found substantially less statistically significant changes in GS end (19% of the stations, mean trend of 1.   Karlsen et al. (2009) found an average trend of − 2.7 days per decade for GS beg , 3.7 days per decade for GS end and 6.4 day per decade for GS len over the Fennoscandia during 1982-2006, of which the GS beg is in good agreement with the other presented estimates of thermally-defined GS beg . Noteworthy, the estimates GS end (and thus GS len ) by Karlsen et al. (2009) are embedded with a considerable amount of uncertainty due to relatively high bias (ranging from − 17 to 14 days) between satellite and field data used to construct the estimates. Therefore, despite the aforementioned studies being fairly well inline with our findings, direct comparison of magnitude of the thermal GS trends between studies is challenging due to differences in analysis periods and methods (i.e. thermal GS vs. satellite remote-sensing), large spatial variability in temporal trends and the coverage of station networks (for thermal GS).

Thermal growing seasons in respect to large-scale processes
The finding that the increase in GS len is mainly due to earlier beginning of GS can be partly explained with the recently established understanding of the changes in atmospheric circulations in northern Europe. According to Räisänen (2019), especially in southern Finland there has been a pronounced negative effect of circulation change for October temperature trends, suggesting change in circulation towards increased frequency of northerly winds. The found smaller temporal trends in GS end in southern Finland and the larger temporal trends in northern Finland are thus in agreement with the result by Räisänen (2019). Moreover, during July-September the effect of circulation change has mainly been positive The y-axes are scaled to unit variance. In the last row, the maps depict the residual variation as a function of latitude (Y, in ETRS89-LAEA) and longitude (X). Effect size (EF) quantifies how much variation a predictor causes on a response variable when analyzed over the station data suggesting more southern air masses to prevail. These results are largely consistent with Nilsen et al. (2017) who found that changes in synoptic circulation have induced warming in July-September and cooling in October in 1981-2010 over northern Europe. However, the observed temperature trends during the growing season can be only partly explained by the changes in the atmospheric circulation (Räisänen 2019). When the effect of circulation is removed, the remaining warming in all months has shown to be robustly positive.
In tandem with longer growing seasons the accumulated GDDS has significantly increased over 93% of the stations. Noteworthy, for many stations the increasing trend e.g. in GDDS has not been linear (see Fig. 1d), but instead an acceleration of the GDDS increase is evident from ca. 1990 onwards. The increase in GDDS has been the most notable in the southwestern part of the domain (Fig. 7d). While the warming throughout the whole GS affects the increase in GDDS, the observed pattern in GDDS trends resonates with the fact that summer temperatures have increased in Europe during the recent decades largely due to anthropogenic forcing (Christidis et al. 2015). In many cases, relatively stable or even reversed GS trends prevailed. Due to this nonlinearity in the observation data our linear trend estimates over the past 70 years are not likely to capture the entire dynamics of recent warming in the region and thus are likely underestimates. The pattern of increase in GDDS near the Baltic Sea also agrees with future projections of the thermal growing season documented in other studies (Ruosteenoja et al. 2011;Zhou et al. 2018). By performing an attribution study, Christidis et al. (2007) detected the human influence on the lengthening of GS. While they used an earlier time period , they also report that the lengthening of GS attributable to human activity is expected to become more prominent in the future as greenhouse gas emissions continue to warm the climate system. In concordance with Christidis et al. (2007) we suggest that anthropogenic global warming has played a significant role in lengthening the growing seasons in northern Europe. A key concern related to increased GDDS and longer growing seasons, particularly GS beg , relates to altered phenology and life cycle of organisms, especially reproduction and voltinism of pest insects in the northern Europe, facilitating to higher outbreak frequency and shifts in distribution (Annila 1969;Bentz et al. 2019).

The role of local environmental conditions in driving GS variability
We found evidence of local environmental mediation on the spatial variation of the GS variables. For example, the used northness predictor (a proxy for incoming direct solar radiation) showed reasonable but relatively weak effects for mean GS conditions-after controlling for the effects of other predictors, growing seasons tend to begin earlier (ca. 1.5 days) and GS being ca. three days longer over southern aspects compared to northern aspects. Slightly stronger relationships were found for forest cover. However, due to the characteristics of our observation data (i.e. weather stations are located on relatively flat ground) the found effects are likely to be underestimates of the true effect, as slope orientation and forest cover have previously been shown to impact spatial variation in various bioclimatic variables (Ashcroft and Gollan 2012;Maclean et al. 2016;Meineri and Hylander 2016;Greiser et al. 2018). This is relevant since for example aspect affects performance of poikilothermic organisms, such as insects (Kantola et al. 2014;Blomqvist et al. 2018). Interestingly, we found relatively strong effects of urban fraction on the investigated GS variables; on average the increase in These results are in agreement with previous research showing how urban climates mediated by land use and built structures can lead to differences in the timing of vegetation phenology compared to surrounding rural areas and within urban areas (Jochner et al. 2012;Zipper et al. 2016). Our spatial modeling further suggests that temporal changes in GS variables have remained smaller close to lakes and sea, and areas of dense forest cover indicating mechanisms related to latent and sensible heat transfer, and radiation interception, respectively, that can have implications on local temperature buffering (Ashcroft et al. 2009;De Frenne et al. 2021). These areas with potentially slower than average velocity of climate change may turn out to be important for e.g. biodiversity conservation under future warming (Loarie et al. 2009;Suggitt et al. 2018;Heikkinen et al. 2020). Similar spatial patterns hold for interannual GS variability (or GS predictability)-our results suggest that there is on average less year-to-year variability in GS variables at high elevation areas, close to water bodies and in areas of high forest cover, although the variability is predicted to increase along urban land use.

Data and modeling uncertainties
Despite the placement of weather stations mostly on lowland, relatively flat and/or open areas (Graae et al. 2012;De Frenne and Verheyen 2016), we were able to produce spatial estimates of the multiple GS variables at unprecedented high-spatial resolution over a large domain. Importantly, these estimates are in good agreement with observations and show meaningful responses to the investigated environmental gradients. Nevertheless, we argue that the inefficiency of weather stations to represent environmental conditions relevant for local ecosystems, often not well represented by the closest weather station or gridded data, is one of the key uncertainties in tracking climate change and its impacts globally. Recently, much effort has been made to develop semi-mechanistic local climate modeling approaches (operating at spatial scales < 0.01 km 2 ) that are based on known physical constraints of air and surface temperature variability (Kearney and Porter 2016;Maclean et al. 2019) and are also applicable over future climate scenarios (Maclean 2020). However, due to computational constraints such models are not yet applicable over large regional domains (such as northern Europe) and over long time periods, and thus empirical approaches built on static landscape variables as proxies for various surface-atmosphere processes are still needed. In addition, other statistical methods for spatial prediction exist ranging from regression to machine-learning (Li and Heap 2011). Geostatistical approaches, such as kriging, have also been widely-used in predicting environmental variables (Matheron 1963;Haylock et al. 2008;Hofstra et al. Fig. 8 The relationships between the spatial variation in temporal trends (β) and the environmental predictors. The fitted functions were estimated using generalized additive models and they depict the relationships between the response variables ( , after holding other predictors constant. The y-axes are scaled to unit variance. In the last row, the maps depict the residual variation as a function of latitude (Y, in ETRS89-LAEA) and longitude (X). Effect size (EF) quantifies how much variation a predictor causes on a response variable when analyzed over the station data 2008). However, in this study the use of e.g. kriging interpolation would not be computationally feasible due to the high-resolution of the analysis. Moreover, Aalto et al. (2013) showed that generalized additive models (GAM, as deployed here) can outperform kriging interpolation in predicting air temperature variability over Finland. Consequently, the large spatial domain poses multitude challenges for the adopted static modeling approach (Meineri and Hylander 2016; Aalto et al. 2017b). For example, the threshold for defining the thermal growing season (here + 5 °C) may not be appropriate for characterizing biological activity in different parts of the study area ranging from Arctic-alpine tundra to southern hemi-boreal forests (Fronzek and Carter 2007). Moreover, sea and lake freeze-up (and melting) dynamics and phenology are likely to substantially differ over the study area. While the estimated relationships between GS and water cover variables are reasonable, such spatially asynchronous seasonal dynamics can increase uncertainty to the estimated relationships. It is also important to note that due to large model extrapolation a fair amount of uncertainty remains in predictions over areas located above ca. 1300 m above sea level. Finally, one deficiency of our approach is that the GS trends were calculated using data with varying temporal coverage. This compromise was necessary in order to increase the station network's environmental coverage which in turn allowed us to predict the spatial variation of the temporal trends over the study domain with reasonable accuracy.

Implications of changing growing seasons and the need for high-resolution GS data
Longer growing seasons and increase in GDDS impact northern ecosystems in various ways, but the issue is still bearing plenty of uncertainty. In general, northern forest ecosystems are resilient to long-term environmental changes, but they are vulnerable to sudden changes due to limitations to adapt promptly to climate driven risks of natural disturbances (Forzieri et al. 2021). Altered tree ecophysiology at high latitudes predispose trees to intensified numbers of biotic agents, like insects and pathogens. Along with forest pests, increasing impact of insect pests in agricultural systems affect production of food, fiber, bioenergy feedstock, and growing renewable materials. According to Lehmann et al. (2020), several current agricultural and forest insect pests are alien invasive species, which benefit from elevating temperatures and contribute to high monetary losses due to mitigation and management. Previously, predicting such impacts has been hampered by the use of climate data having spatial resolution not relevant for the modeling targets (Potter et al. 2013). A benefit following longer growing seasons and a warming trend is evidently seen in higher production and spatiotemporal shifts in cultivation areas of agricultural crops (Peltonen-Sainio and Jauhiainen 2020) and in stem volume increment of trees (Kauppi et al. 2014;Henttonen et al. 2017). However, the impact of climate-driven natural disturbance agents increasingly cut crop and forest biomass production and lower the profit. The carbon sink of forests may not persist continuously in the long term (Nabuurs et al. 2013). Above-mentioned cases emphasize the need of modeling tools targeted for accurate predictions in agricultural and forest ecosystems, applying high-resolution gridded data as a source for modeling approaches. Detailed thermal growing season data at 100 m resolution provide valuable information, e.g., for ecosystem and risk models (Netherer et al. 2019), forest growth (Minunno et al. 2019), and insect migration models (Zurell et al. 2016). Ultimately, such data provide us with means for mitigation and adaptation to changing climate, and supervising decisionmaking and management options. For example, knowledge on growing seasons is vital in agriculture (Peltonen-Sainio et al. 2016), horticulture (Campoy et al. 2011) and forest management (Castanha et al. 2013), both in theory and in practice. In addition, plant breeding uses detailed information on variation of growing seasons (Lillemo et al. 2010) and practitioners are well-familiar with the concept of GDDS. Climate warming affects natural disturbances of forest trees such as the life cycle and voltinism of bark beetles (Bentz et al. 2019).

Conclusions
We found remarkable spatio-temporal variability in GS variables over northern Europe. Such variation was strongly controlled by latitudinal and elevational gradients albeit partly constrained by local topography, proximity to water bodies, forest cover and urban land use. By deploying a large network of weather stations this study provides the most up-to-date and spatially detailed understanding of the current variation and past changes in multiple GS variables in northern Europe. Past changes in GS variables could be used as predictors for temporal and spatial shifts in occurrence of various organisms. As thermal growing season conditions are the most proximally related to the organisms' performance, the produced high-resolution data layers will be highly relevant in various ecosystem applications such as in assessing suitable conditions for crop production and modeling outbreak risks for forest insect pests.

Data availability
The main output data layers (i.e. mean GS conditions and temporal trends) are deposited in the public repository of the Finnish Meteorological Institute at https:// doi. org/ 10. 23728/ fmi-b2sha re. ff569 ebaa6 2a4e7 5baca fd76a d2a9a b6.
Code availability The R code used for producing the results is available from the corresponding author by request.

Conflict of interest
The authors declare no conflict of interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.