Evapotranspiration Intensification Over Unchanged Temperate Vegetation in the Baltic Countries Is Being Driven by Climate Shifts

Shifts in climate driven by anthropogenic land use and land cover change are expected to alter various land–atmosphere interactions. Evapotranspiration (ET) is one of these processes and plays a fundamental role in the hydrologic cycle. Using gridded reanalysis and remote sensing data, we investigated the spatiotemporal trends of precipitation, temperature, and ET for areas in the Baltic countries Lithuania, Latvia and Estonia where the land cover type had not changed from 2000 to 2018. We focused on ET but investigated the spatiotemporal trends for the three variables at monthly, seasonal, and annual time scales during this period to quantify trade-offs among months and seasons. We used the Mann-Kendall test and Sen’s slope to calculate the trends and rate of change for the three variables. Although precipitation showed fewer statistically significant increasing and decreasing trends due to its high variability, temperature showed only increasing trends. The trends were concentrated in late spring (May, +0.14°C annually), summer (June and August, +0.10°C), and early autumn (September, +0.13°C). For unchanged forest and cropland areas, we found no statistically significant ET trends. However, Sen’s slope indicated increasing ET in April, May, June, and September for forest areas and in May and June for cropland. Our results indicate that during the study period, the temperature changes may have lengthened the growing season, which affected the ET patterns of forest and cropland areas. The results also provide important insights into the regional water balance and complement the findings of other studies.


INTRODUCTION
Changes in precipitation and temperature regimes differ among the world's climatic regions. Previous studies have demonstrated changes or high variability in the patterns of two key climate variables, precipitation and temperature, in recent years in a range of climatic regions (Van den Silva Junior et al., 2018;Haghtalab et al., 2019;Krauskopf and Huth, 2020). The impact of these changes, which varies among the climatic regions, is exemplified by an increased number of extreme events, such as heat waves, intense rainfall, or long-term droughts (Sena et al., 2012). In addition, according to the Intergovernmental Panel on Climate Change (IPCC) (Seneviratne et al., 2012), the frequency and magnitude of these extreme events and the climate variability are expected to increase in the future.
Among the factors responsible for the changing climate, anthropogenic land use and land cover (LULC) change is a main factor (Pielke, 2005;Ceccherini et al., 2020;Sy and Quesada, 2020), and is having a stronger impact than was assumed in previous studies (Seneviratne et al., 2012). Depending on the extent and intensity of the LULC change, climate variables can be affected at local, regional, and global scales (Hale et al., 2006;Mahmood et al., 2014;Cao et al., 2015;Chen and Dirmeyer, 2017;Findell et al., 2017). In Eastern India, for instance, 25 to 50% of the overall temperature increase of ∼0.3 • C from 1981 to 2010 could be explained by local LULC change (Gogoi et al., 2019). In the Amazon, Llopart et al. (2018) analyzed how the spread of deforestation affected the regional climate using simulated scenarios and climate models. They reported that deforestation caused a dipole pattern in response to the region's precipitation, with increased precipitation in the eastern Amazon (+8.3%) and decreased precipitation in the west (−7.9%). At a global scale, Lawrence and Chase (2010) found that anthropogenic land cover change imposed a warming trend in the near surface atmosphere and reducing precipitation regionally, and both climate factors varied greatly between regions of the world.
In most cases, the impact of LULC change on climate results from changes in the biophysical and biogeochemical processes that regulate the complex land-atmosphere interactions (Pongratz et al., 2010;Alkama and Cescatti, 2016;Perugini et al., 2017). Evapotranspiration (ET) is one of the processes affected by climate and LULC change, and it plays an important role in balancing the hydrological cycle (Zhang et al., 2016). In temperate forest regions, for example, ET returns 71% of the precipitation to the atmosphere (Miralles et al., 2011). Using ET data derived from remote sensing or modeling, studies have been exploring the consequences of climate and LULC change on the dynamics of ET (Dias et al., 2015;Spera et al., 2016;Poon and Kinoshita, 2018;Paca et al., 2019;Qu and Zhuang, 2019;dos Santos et al., 2020). Wang et al. (2019) analyzed the spatiotemporal variability of ET over a Chinese catchment that has experienced strong climate and LULC changes. They found that the annual mean ET increased by nearly 14 mm as the result of a warming and drying trend combined with increasing areas of forest, grassland, water bodies, and urban land. Using a water-balance approach (i.e., the difference between precipitation and stream discharge) to estimate ET, Hamilton et al. (2018) found that ET is resilient in temperate humid catchments that have been experiencing climate and LULC change; that is, ET remained relatively stable. They attributed this resilience to field measurements that showed the ET rates did not differ greatly between rain-fed annual crops and perennial vegetation, and to the fact that these LULC types showed the biggest changes in area during the study period. However, this approach was limited because the ET was only calculated for the entire catchment and not independently for the different LULC types.
Even though several studies have provided spatiotemporal analysis of the combined effects of climate and LULC changes on ET (Mao et al., 2015;Li et al., 2017;Teuling et al., 2019), fewer studies have investigated the changes in ET in areas where LULC has not changed (Gaertner et al., 2019), where changes in ET are most likely driven by changes in climate variables. The limited number of studies that have investigated ET changes in areas where LULC has not changed, such as many forested areas, have reported different patterns for different climate regions. For semi-deciduous tropical forests, Vourlitis et al. (2014) showed that ET has decreased, whereas for temperate forests, Gaertner et al., 2019) reported that ET increased due to the longer growing season. Most of these studies evaluated annual ET changes, and restricted the assessment of the trade-offs between monthly and seasonal time scales (e.g., increasing ET in spring can be offset by an ET decrease in summer). Moreover, the studies did not account simultaneously for spatial variation and changes over time. Therefore, they did not produce a clear understanding on how ET changes spatially and temporally at different time scales (e.g., monthly, seasonally, and annually) in areas with stable LULC. Studies of ET in this context of stability are important because they help us to understand how stable LULC types such as forest and agricultural areas have been adapting to the new climate conditions and reveal the importance of these areas as a long-term solution to minimize climate change by stabilizing water cycles and sequestering carbon (Funk et al., 2019).
Given the abovementioned limitations of existing research and the lack of understanding of how ET rates change spatially and temporally at large scales, we designed the present study to investigate the changes in ET over areas with stable forest and cropland areas. We focused our analysis on three Baltic countries (Estonia, Latvia, and Lithuania). We used a newly developed gap-filled ET remote sensing dataset (Running et al., 2019) that provided data from 2000 to 2018. We chose these countries because the region has seen a clear shift in its climate from 1951 to 2015, with gradual warming and a slight trend of increasing precipitation (Jaagus et al., 2017). In addition, most of the region has not shown significant LULC changes during the study period, and no similar study has been performed in these three countries.
We defined unchanged forest as areas that showed no significant disturbance caused by anthropogenic activities such as clear-cut or selective logging during the study period. We defined unchanged cropland as cultivated land that remained cultivated throughout the study period but did not consider the effects of crop changes (e.g., from winter wheat to corn) to represent an LULC change. We chose forest and cropland areas for our analysis because these LULC types cover the majority of Estonia, Latvia, and Lithuania, and both types show large continuous patches of unchanged LULC. We also explored the spatiotemporal patterns of the two key factors that affect ET (i.e., precipitation and temperature) using gridded reanalysis data to better understand the driving forces for the changes in ET.
We posed two questions: (i) What are the spatial and temporal patterns of temperature, precipitation, and ET throughout the Baltic countries? (ii) How has ET of unchanged forest and cropland areas changed in this region? We hypothesized that even though the region has not experienced intense LULC change since 2000, the overall ET rates of forest and croplands have changed to adapt to the observed shifts in the climate regime in response to large-scale atmospheric circulation patterns (Jaagus et al., 2003). Note that although climate change is the terminology adopted by most similar studies, we decided to use the term climate shift because our study period is shorter than the period (at least 30 years) recommended for climate change studies (Seneviratne et al., 2012).

Study Area
The study area encompasses 175 × 10 3 km 2 of the three Baltic countries in north-eastern Europe. The region's climate is characterized as humid continental Dfb according to the Köppen climate classification, with mean annual temperature ranging between 5.5 and 6.5 • C and mean annual precipitation ranging between 600 and 750 mm. In the westernmost part of the study region, which lies on the coast of the Baltic Sea, winters are much milder and the Köppen climate type changes to Cfb. However, parts of the region have shown shifts in the precipitation and temperature trends. The precipitation has increased during the cold months of the year (November to March) and in June, but the annual mean temperature has increased by 0.3 to 0.4 K per decade from 1951 to 2015 (Jaagus et al., 2017).

Precipitation and Temperature Data
We used the gridded daily precipitation and daily mean temperature products from version v20.0e of the European Climate Assessment and Dataset, which is called E-OBS (Cornes et al., 2018). E-OBS was developed by interpolating climate data between weather stations based on data provided by the national meteorological agencies of the European countries. The product is available as a regular grid with a resolution of 0.1 or 0.25 • , and covers the period since 1950. We used the 0.1 • -resolution dataset, which represents a resolution of about 10 km. E-OBS has been applied in climate monitoring studies (Nastos et al., 2013;Van Der Schrier et al., 2013) and for validation of climate products (Katsanos et al., 2016;Domínguez-Castro et al., 2020). The reported shortcoming of this dataset is the uneven distribution of weather stations, which compromises the interpolation of E-OBS precipitation and temperature estimates (Hofstra et al., 2010). However, this limitation is less problematic in our study region, which has a higher weather station density (Navarro et al., 2019) than in other areas of Europe.
Even though the E-OBS grid was developed from national quality-controlled data, we checked for potential errors (e.g., missing or extreme low and high values) in the gridded precipitation and temperature datasets. We found no errors. To perform our analysis, we aggregated the gridded 0.1 • daily precipitation and temperature data into monthly, seasonal, and annual temporal scales. We defined the seasons using the standard periods used for northern hemisphere climatology (Jaagus et al., 2017), with spring comprising March, April, and May; summer, June, July and August; autumn, September, October, and November; and winter, December, January, and February. For precipitation, we summed all total daily values for each month, season, and year, but for temperature, we calculated the mean for each temporal scale.

Evapotranspiration Data
To perform the ET analysis, we used the relatively recently developed MODIS MOD16A2GF product that comprises gapfilled 8-day composite images produced at 0.00416 • pixel resolution (∼500 m) (Running et al., 2019). The gap-filled data is generated by mainly removing contaminated FPAR/LAI input data. For the cases where the pixel of FPAR/LAI pixels did not meet the quality criteria, the value is determined by linear interpolation between the previous period's value (Running et al., 2019). The MODIS ET pixels represent the sum of ET (kg/m 2 ) for 8-day periods. The product is an improvement of the MOD16A2 product that eliminates poor-quality pixels (e.g., pixels contaminated by cloud cover or aerosols). The MOD16A2 product provides almost real time estimates of ET using the Penman-Monteith equation, and has been used in several studies (Spera et al., 2016;Moura et al., 2019). We acquired 8-day composite MOD16A2GF images from January 2000 until December 2018.
Before analysis, we removed MODIS ET pixels with values ranging from 32,761 to 32,767 by assigning those values of NoData (the coding used by the R and QGis software). These pixels are flagged in the original dataset as non-vegetated areas, and therefore, no ET rate was calculated (Running et al., 2019). Moreover, we removed all pixels that showed negative values. In this way, we only used pixels that had plausible values during the entire period of analysis (2000 to 2018). Because of the pixelremoval process, our ET analysis covered 70% of the total Baltic countries territory (∼120 × 10 3 km 2 ). To retrieve the monthly ET rates, we calculated the mean of the 8-day composites images within the respective month for each pixel. The mean value was then divided by eight to obtain daily values, which were finally multiplied by the number of days of each month. We used this approach to retrieve the monthly ET rates because April, May, October, and November had different number of 8-day composite images in leap years. The seasonal and annual ET rates were calculated by aggregating the monthly values. We did not analyze the winter season for four reasons: (i) No data was available in January and February of 2016. (ii) The winter months of other years showed many pixels with negative values, and this would have resulted in the exclusion of a large number of pixels since we worked only with pixels that had plausible values for the whole period. (iii) ET rates in winter are low and the values are negligible, therefore studies have mainly used only the months of the growing season (Gaertner et al., 2019;Niu et al., 2019;Ruiz-Pérez and Vico, 2020). (iv) The ET analysis was not stratified by forest type and as deciduous and evergreen forest are present in the region, the ET rates analysis would be bias since some forest types loses or maintain the leaves. Although the original unit for the ET rates was kg/m 2 , we decided to convert this to mm (at a rate of 1 kg/m 2 for 1 mm), to align with the units used in previous studies.

Land Use and Land Cover Data
For LULC data, we used version 2.0.7 of the (CCI Land Cover product 1 ). This product consists of annual time series of consistent global land cover maps with 300-m spatial resolution from 1992 to 2018. We used the product data from 2000, 2005, 2010, 2015, and 2018. For 2018, we used version 2.1.1 of the product, which is consistent with v2.0.7. The advantage of the newer dataset is its consistency over time. We reclassified the original 22 LULC types into 6 main classes based on the IPCC class definitions (ESA, 2017 Land Cover CCI Product user guide version 2.0, 2017) (Supplementary Table 1). The final LULC types were cropland, forest, grassland, urban, wetland, and other. We identified the pixels of forest and cropland that remained as these classes throughout the study period to extract the unchanged forest and cropland areas. The map was only created for forest and cropland because these two LULC types covered more than 87% of the entire study region. Moreover, these LULC types also showed large continuous areas in which the LULC did not change, unlike in the other LULC types, and this lets us test our research hypothesis.

Analysis
First, we examined the monthly, seasonal, and annual spatial and temporal trends for temperature and precipitation for the entire region using the E-OBS data, and the MOD16A2GF data for ET. We examined the trends at a pixel level, with a spatial resolution of 0.00416 • . Second, we assessed the monthly, seasonal, and annual ET for unchanged forest and cropland areas using the MOD16A2GF data. We used the Mann-Kendall test to detect significant increasing or decreasing trends.

Spatial Analysis
To analyze ET over forested and cropland areas, we calculated the vegetation cover of forest and cropland within the MODIS ET pixels that were overlapped by these LULC types. The MODIS ET pixels were defined as forest or cropland pixels if, during the entire study period, these LULC types covered >50% of the pixel area. The MODIS pixels that showed <50% of its area as forest or cropland during the study period were considered to have changed, and were excluded from this analysis. In this way, we could investigate the MODIS ET pixels that showed a decrease, increase, or no trend during the study period for areas that had forest or cropland as the major (>50%) LULC type.
To obtain the most representative ET values for unchanged forested and cropland areas, we created a buffer that extended 1 km inward from the edge of the unchanged area and extracted only those pixels that were completely within the unchanged area. This removed the potential bias due to the edge effect for the border pixels. After these procedures, we retrieved 15,359 pixels of forest and 899 pixels of croplands as the most representative MODIS ET pixels of these LULC types. We averaged the monthly, seasonal and annual ET values of these pixels for the trend analysis. 1 https://www.esa-landcover-cci.org/?q=node/164

Trend Analysis
We used the non-parametric Mann-Kendall test (Mann, 1945;Kendall, 1957;Khambhammettu, 2005) to identify pixels with significant temporal trends (significant increase or decrease, or no significant trend) for the three variables (precipitation, temperature, and ET). The Mann-Kendall test has been widely applied for trend analysis of time series of meteorological (Jaagus, 2006;Almeida et al., 2017;Atta-ur and Dawood, 2017) and hydrological data (Hamed, 2008;Dinpashoh et al., 2011;Jaagus et al., 2017). However, the Mann-Kendall test has been shown to be sensitive to autocorrelation in the input data (Kumar et al., 2009). Therefore, before applying the Mann-Kendall test, we calculated the lag 1 coefficient of correlation test for serial correlation (Bari et al., 2016;Gao et al., 2020) in the time series of precipitation, temperature and ET for all time scales of the analysis. The calculation indicated that an average of <3% of all pixels analyzed showed significant autocorrelation in all-time series. Thus, we decided not to apply any correction and to use the Mann-Kendall test on the original data. We used a confidence level of 95% to define significant decreases and increases in the time series for each pixel (p < 0.05).
To estimate the magnitude of the trend, we used the nonparametric Sen's slope value (SS) method (Sen, 1968). The magnitude of the change is given by the median slope value for all pairwise points in the time series. Negative values of SS indicate a decreasing trend, whereas positive slopes indicate an increase. To calculate the Mann-Kendall trend and SSs, we used the package "wql" (Jassby and Cloern, 2017) implemented in version 3.6.1 of the R software. We calculated the trends as well as the SS values for all pixels at monthly, seasonal, and annual scales. We calculated SS values even for those pixels that did not show significant trends according to MK test.

Land Use and Land Cover Change
In 2000, 53.1% of the Baltic countries were covered by forest, followed by cropland (36.7%) and grassland (5.2%) (Supplementary Table 2). Urban, wetland and other land together totaled 5.0% of the region (Supplementary Table 2). By 2018, forest area had decreased by 0.4% and cropland had increased by 0.4% (Supplementary Table 2). Although the total change in the area of LULC types was negligible, the changes of areas of cropland and forest to other classes were noticeable (Figure 1). For example, 2,532.1 km 2 of forest was converted into cropland and 801 km 2 into grassland from 2000 to 2018. At the same time 2,300.8 km 2 of cropland was converted into forest (Figure 1). Despite these changes, 72.8% of the forest area and 73.2% of the cropland area in 2000 remained completely unchanged throughout the study period (Supplementary Table 3).

Spatiotemporal Precipitation Analysis
Our spatial and temporal analyses of the monthly precipitation at a pixel level indicated that half of the months showed no significant trend in the 2,557 pixels analyzed between 2000 and 2018 (Figure 2A). April, September, and December showed a significant increase in 13.6, 3.1, and 8.1% of the total pixels (Supplementary Table 4), respectively. March, June, and July presented significant decreasing trends in a small proportion of the pixels (0.8, 1.4, and 0.2% of the total, respectively). Spatially, it is interesting to note that the pixels with the significant decreasing trend were mainly concentrated in the south-west of Lithuania (in March and June), whereas the increase in precipitation was not only in the south and southwest of Lithuania and Latvia but also in the northeast of Estonia (Figure 2A). The SS values (Figure 2B) revealed the magnitude of the change for the pixels with a significant trend (p < 0.05) and without a significant trend (p > 0.05). September showed an average SS value of 2.37 ± 0.29 mm per year (mean ± SD) for the pixels with a significant increase, for a total increase of 45 mm per pixel for the whole study period. The pixels with an increasing trend in April and December showed an average SS increase of 1.60 ± 0.26 and 1.31 ± 0.32 mm per year, respectively (Supplementary Table 4). Similar SS values were observed for the pixels with a significant decreasing trend in March, June, and July, with average SS decreases of 1.66, 2.33, and 2.22 mm per year, respectively (Supplementary Table 4).
Although the trend analysis for seasonal precipitation (Supplementary Figure 1) showed that spring had only four pixels with a significant increasing trend and 31 pixels with a significant decreasing trend, the average SS reveals an overall increase trend of 0.32 ± 1 mm per year (mean ± SD) in the region. The pixels with an increasing trend in spring had an average SS value of 1.40 ± 0.24 mm per year, versus −1.34 mm per year for pixels with a decreasing trend (Supplementary Table 5). Winter pixels with an increasing trend (7 piels) had a higher average SS value (5.63 mm per year) than the spring pixels with an increasing trend (1.4 mm per year). Autumn and summer showed no pixel with significant trends. Despite the lack of significant trends, summer was the only season that showed a decrease in precipitation for the pixels that showed no significant trend (mean SS = −0.91 mm per year). The annual analysis revealed no significant trend for any pixel (Supplementary Figure 2 and Supplementary Table 6).

Spatial and Temporal Temperature Analysis
Analysis of the monthly mean temperature showed that 8 months presented no significant (p < 0.05) trend in any pixel ( Figure 3A). Moreover, no pixel with a significant decreasing trend was observed in any month. Pixels with significant increasing trends were identified in May (6.2% of the 2,557 pixels), June (4.6%), August (16.6%) and September (89.3%). In September, the pixels with a significant increasing trend covered almost the entire region. However, in May, June and August, the pixels with an increasing trend were clustered into a few areas in northern and western Estonia and eastern Latvia and Lithuania. The average SS values of the pixels with a significant increasing trend ( Figure 3B) were 0.14 ± 0.14 • C per year in May (mean ± SD), 0.10 ± 0.01 • C per year in June and August, and 0.13 ± 0.02 • C per year in September (Supplementary Table 7).
In the seasonal analysis, only autumn had pixels with significant increasing trends (6.6% of the pixels), and these pixels were located in the west of Lithuania and Latvia (Supplementary  Figure 4). The pixels with an increasing trend in autumn had SS values that showed an average increase of 0.11 ± 0.01 • C per year (Supplementary Table 8). Among the pixels with no significant trend in the other seasons, winter presented the highest average SS (0.07 • C per year) with values ranging from 0.14 to −0.02 • C per year (Supplementary Table 8). The annual analysis showed that 31.8% (813 out of 2557) of the pixels presented increasing trends (Supplementary Figure 5) and they are clustered in the west and east regions of Latvia and Lithuania, and in a smaller area in northern Estonia. The average SS of the pixels with an increase was 0.07 • C per year, with values ranging from 0.04 to 0.09 • C per year, versus 0.05 • C per year for pixels with no trend and values ranging from 0.02 to 0.08 • C per year (Supplementary Table 9).

Spatial and Temporal ET Analysis From 2000 to 2018
Pixels with significant (p < 0.05) increasing and decreasing trends for ET ( Figure 4A) were more scattered compared to pixels with significant temperature and precipitation trends ( Figure 3A). In May,June,and September,23.2,9.1,and 15.2% of the pixels (from a total of from 1 423 312) showed increasing ET trends, with mean SS values of 1.41 ± 0.50 (mean ± SD), 1.69 ± 0.67, and 0.65 ± 0.31 mm per year (Figure 4B), respectively. Interestingly, the mean SS values for the pixels with a decreasing trend in May, June, and September were much less common (0.1, 0.2, and 0.1%, respectively), showed SS values similar to those of the pixels with an increasing trend (Supplementary Table 10). August showed spatial segregation, with pixels that showed an increasing trend (9.3% of pixels) in the eastern part of the region and decreasing trends (3%) in the western part of the region ( Figure 4A). November showed the biggest percentage of pixels with a decreasing trend for ET (9.5% of the total pixels; Supplementary Table 10). However, the mean SS values were smaller than the values for pixels with an increasing trend (Supplementary Table 10). March, April, July, and October showed <2.5% of the pixels with either increasing or decreasing trend. For these months (except July), the pixels with no significant trend revealed an average increase in the ET rates (Supplementary Table 10).
By analysing the ET trends seasonally, we found that 13% of the pixels in autumn (Supplementary Table 11) and 13.5% in spring presented a significant increasing trend, which is more than the proportion in summer (10%). All seasons had fewer than 1% of the pixels showing a decreasing trend (Supplementary Table 11). The pixels with a decreasing trend were spread throughout the Baltic countries during summer and autumn (Supplementary Figure 6). We also observed the latter spatial pattern for pixels with an increasing trend in all three seasons. On an annual basis, the spatiotemporal analysis revealed that 34.9% of the pixels had a significant increasing trend and only 0.3% showed a significant decreasing trend (Supplementary Figure 7 and Supplementary Table 12).
We found that 83.5% of the MODIS ET pixels had forest (48.1%) and cropland (35.4%) as the major LULC type (with these types accounting for >50% of the MODIS ET pixel area). Hereafter, we refer to these as forest and cropland pixels. Most of these pixels (>77%) had no significant trend in ET in any month. Few forest and cropland pixels (<10%) showed a significant decreasing trend, with August showing the highest percentage (5.76%) for forest pixels (Figure 4C). For pixels with increasing trends, pixels of both LULC types showed similar percentages for March, April, October, and November (<2%). However, the difference between crop and forest pixels with an increasing trend was <4 percentage points in May, June, July, and August, September had a difference of 7.3 percentage points (11.2% of forest pixels and 3.8% of crop pixels; Figure 4C).

Evapotranspiration Over Unchanged Forest and Cropland Areas
We extracted the monthly ET of 10 949 MODIS ET pixels dominated by unchanged forest and 377 pixels dominated by unchanged cropland. Using the average monthly ET of these pixels, the Mann-Kendall trend analysis revealed that neither forest nor cropland showed a significant increasing or decreasing  for cropland also showed increase for most months, except August (−0.2 mm per year) and November (−0.01 mm per year). Crop ET in June presented the highest value (0.64 mm per year), followed by May (0.50 mm per year). On a seasonal basis, forest showed a significant increasing trend in autumn (SS = 0.46 mm per year). Moreover, forest showed higher SS values than cropland in all seasons (Supplementary Table 14). At an annual scale, only forest had significant increasing trend (3.46 mm per year).

Precipitation and Temperature
In general, our results agreed with those of previous studies (Kriauciuniene et al., 2012;Jaagus et al., 2014Jaagus et al., , 2017, which also identified shifts in precipitation and temperature at monthly, seasonal, and annual temporal scales using weather station data. However, our study advanced these studies of climate shifts by analysing them both spatially (i.e., using a gridded dataset for the whole region) and at three temporal scales. The monthly analysis provided new insights that analysis with a coarser temporal resolution (i.e., seasonal and annual) would have not provided. For instance, the SS values for the overall summer season showed a reduction in precipitation. However, the monthly analysis showed that June and July had areas with statistically significant decreasing trends for precipitation and negative mean SS values for the pixels with non-significant trends, whereas August showed positive SS values. Therefore, the negative SS values of the summer season were influenced mostly by decreased June and July precipitation but not by August changes. This highlights the necessity of evaluating trends at a monthly level, as this more clearly reveals differences between the months within a season. Moreover, the spatial analysis using the monthly time scale provided an overview of the reduction in precipitation, which was more spatially dispersed in June and more clustered in July and August. Similar observations can be done for spring. March and May (the beginning and end of spring) showed reduced precipitation (negative SS values), whereas April (in the middle of the spring) showed an increase in precipitation, although with no significant trend. This dispersion of statistically significant trends at monthly and seasonal level was also reported by Jaagus et al. (2017), who investigated precipitation trends in Estonia using weather station data. They suggested that largescale atmospheric circulation patterns were the main factor for the high spatiotemporal variability in the precipitation over the region (Supplementary Figure 3).
As was the case in other temperate regions Oishi et al., 2018), our results also indicated that the Baltic countries are experiencing a warming shift in most months, in all seasons (although only autumn showed a significant increasing trend), and in the mean annual temperature. Most of the statistically significant increase in temperature at a monthly scale was detected at the end of spring (May), at the beginning (June) and end of summer (August), and at the beginning of autumn (September). Similar findings were reported by Kriauciuniene et al. (2012), although these authors analyzed a longer time series record (from 1922 to 2007) for Lithuania, Latvia and Estonia.
By combining the temperature and the precipitation analysis, we found critical periods. For example, at the end of spring (May) and beginning of summer (June), the temperature showed statistically significant increasing trends, whereas the precipitation during the same period showed statistically significant decreasing trends. This can limit soil water availability during the spring warming period, and this will have several consequences for plant development, especially during the growing season. According to Grillakis (2019), soil moisture droughts events are expected to increase in frequency and spatial extent, and eastern Europe is likely to be most strongly affected. Similarly, Samaniego et al. (2018) reported that up to 40% of European soils will be affected by soil drought caused by anthropogenic warming if the temperature increases by 3 • C compared to the pre-industrial period.

Evapotranspiration
Based on the ET trends and the LULC change analysis, we confirmed that the changes in ET rates over the Baltic countries are most likely driven more by shifts in the climatic variables than by LULC change. This assumption agrees with the results of Li et al. (2017), who demonstrated a more significant impact of climate change than LULC change on ET throughout China. These changes in ET have an impact on the regional water balance. As highlighted by Teuling et al. (2019), the reduction in precipitation combined with an increase in ET to reduce streamflow. For example, by analysing precipitation, temperature and river discharge trends from 1991 to 2007, Kriauciuniene et al. (2012) found a significant decrease in river discharge in the western part of the Baltic region during the spring, a decrease in the whole region during the autumn, and an increase throughout the region during winter. Our results for ET showed statistically significant increases for most months that showed a reduction in river discharge in the same study area, and this partially explains the decrease in the river discharge. However, the reduction in the number of days and depth of snow cover (Rimkus et al., 2018) between 1961 and 2015 also affected the regional water balance, especially by limiting the amount of snowmelt water and groundwater, which are the major water sources for several rivers (Kriauciuniene et al., 2012). Moreover, our monthly ET analysis indicated that May, June, and July are critical months, since ET exceeds precipitation, on average, in 77% (+31.1 mm), 74.3% (+52.5 mm), and 50.7% (+36.3 mm), respectively (Supplementary Figure 8). This can potentially result in a reduction of soil moisture and groundwater availability if the water transported into the atmosphere by ET is not returning to the soil as precipitation, especially in June and July, when precipitation was decreasing. This type of seasonal hydrological shift was reported as a factor for summer drought in temperate forests (Parida and Buermann, 2014).
The changes in ET for croplands and forested areas are potentially a consequence of the lengthening growing season in the study region. Gaertner et al. (2019) found that an increase of ∼22 days in the growing season since 1982 until 2012 resulted in an increase of 12 mm in ET in the temperate forest region. They reported that temperature, vapor pressure deficit, wind, and humidity are important climatic variable for the change of the growing season. In addition, precipitation has been reported as an important factor for increasing growing season length (Dragoni and Rahman, 2012), since higher growingseason precipitation allows the vegetation to remain active over a longer period (Hussain et al., 2019). Therefore, the trends of increasing precipitation in April, August, and September (significant and not significant) combined with the temperature trends in May, August, and September that we observed, have likely driven the longer growing season. A longer growing season has been already reported for the study area using analysis of daily mean air temperature (Linderholm and Walther, 2008). According to the authors, the growing season increased between 11.5 and 12.2 days mainly due to an earlier starting of ∼10 days and a later ending of ∼2 days. Considering the findings of Gaertner et al. (2019), that a longer growing season result in an increase of ET, we can suggest that the forested ET rate increase in April and September (Figure 5) and the croplands ET rate increase in May (Figure 5) are likely caused by the changes in the growing season duration over the study area. We found that May and June showed reductions in precipitation and increases in temperature. This combination can potentially affect crop development, resulting in a soil moisture decrease that can cause drought during the early stages of crop development, and this can potentially result in significant yield reductions. For instance, the combination of decreased precipitation with decreased soil moisture can reduce crop yield (for wheat, canola, and barley) by 25 to 45% in a dry growing season (Madadgar et al., 2017).
For the forested areas, on the other hand, the reduction of ET during June and July, although not statistically significant, can be related to the decreased precipitation during these months. As a consequence, the forested areas can become more vulnerable to fires (Donis et al., 2017). One of the reasons that the decreases in ET over forested areas in June and July were not significant, and were lower than the decrease over the cropland, is that the forest can recycle more moisture due to its deeper roots (Zha et al., 2010). Lian et al. (2020) found that spring and summer soil moisture in the northern hemisphere has decreased from 1982 to 2011 due to earlier vegetation greening. They observed that this process was driven by an increase in ET rather than increased precipitation. Our observations for forest and cropland areas also agree with the results of Rimkus et al. (2017), who used the normalizeddifference vegetation index, the vegetation condition index, and the standard precipitation index to study drought in the Baltic countries, concluded that cropland was affected more strongly than forest by precipitation deficits.
There are certain limitations of our study. First, our observations for unchanged cropland areas are limited by the small size of the patches of this LULC class, which did not allow us to sample <1.2% of MODIS ET pixels. To better estimate the changes in such small patches, future studies would benefit from higher resolution ET and climate datasets (Ruiz-Pérez and Vico, 2020). The datasets could be created using modeling (e.g., for ET, using the Penman-Monteith or Priestley-Taylor models). By using a higher spatial resolution, we could expand our analysis to include other LULC classes. For example, although we consolidated all forest types into a single forest class and all cropland types into a single cropland class, forests with different species or degrees of vegetation cover and different crops have different ET characteristics. In future research, these differences should be accounted for in our analysis. Another shortcoming is the lack of an uncertainty assessment for the MODIS ET data over the study region. However, several studies have pointed out the reliability of this data (Aguilar et al., 2018;Sriwongsitanon et al., 2020), with very good results over forested areas (Paca et al., 2019) and cropland areas (Velpuri et al., 2013). Future studies should also investigate the relationship or combined effect of precipitation and temperature change on ET rates at different time scales. Furthermore, it would be interesting to utilize water balance models to investigate how the changes in ET are affecting the water balance of the region (Kumar and Merwade, 2011).

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.