Soil Moisture Variations From Boreal Forests to the Tundra

Soil moisture has a profound influence on life on Earth, and this vital water resource varies across space and time. Here, we explored soil moisture variations in boreal forest and tundra environments, where comprehensive soil moisture datasets are scarce. We installed soil moisture sensors up to 14 cm depth at 503 measurement sites within seven study areas across northern Europe. We recorded 6,138,528 measurements to capture soil moisture variations of the snowless season from April to September 2020. We described the spatio‐temporal patterns of soil moisture, and test how these patterns are linked to topography and how these links vary in space and time. We found large spatial variation and often less pronounced temporal variation in soil moisture across the measurement sites within all study areas. We found that throughout the measurement periods both univariate and multivariate models with topographic predictors showed great temporal variation in their performance and in the relative influence of the predictors within and across the areas. We found that the soil moisture‐topography relationships are site‐specific, as the topography‐based models often performed poorly when transferred from one area to another. There was no general solution that would work well in all the study areas when modeling soil moisture variation with topography. This should be carefully considered before applying topographic proxies for soil moisture. Overall, these data and results highlight the strong spatio‐temporal heterogeneity of soil moisture conditions in boreal forest and tundra environments.

Topography is a static factor. Thus, relating it directly to such a dynamic phenomenon as soil moisture can be problematic, and yet this is often done Riihimäki et al., 2021), for instance, in ecological and biogeographical modeling (see e.g., Niittynen et al., 2018). Topography-based variables, such as the Topographic wetness index (TWI), are commonly used proxies for soil moisture (Ågren et al., 2014;Kopecký et al., 2021;Sorensen & Seibert, 2007;Western et al., 1999). Topography-based indices show temporal variation in their performance as proxies for soil moisture; they perform better in wetter seasons (Ali et al., 2014;Riihimäki et al., 2021;Western et al., 1999). Topography is most strongly related to runoff (Beven & Kirkby, 1979), and to some extent, topography also influences the local spatial variation in infiltration (via e.g., its effects on soil formation), evaporation (via microclimate), and transpiration (via vegetation). All these processes vary greatly also in time. However, topography-based indices do not contain any temporal information relevant for soil moisture, for example, snow melt, precipitation, evaporation, or transpiration (Crave & Gascuel-Odoux, 1997;Western et al., 1999;Wilson et al., 2004). Overall, it is insufficiently understood how the soil moisture-topography relationship varies in different environmental contexts and points in time.
Here, we explore surface soil moisture across and within boreal forest and tundra environments. We use temporally continuous soil moisture sensor data (up to 14 cm depth) from 503 measurement sites within seven study areas. The data range from April to September, which is the snowless season in the study areas. The objectives of this work are to (a) describe the spatio-temporal patterns of soil moisture across and within boreal forest and tundra environments, and (b) test how these patterns are linked to topography and how these links vary in space and time. Our aim is to describe the variability of soil moisture and its relationship with topography across time and space.

Study Areas
The seven study areas extended from southern Finland to northern Norway (Table 1, Figure 1). The areas covered distinct Fennoscandian environments, that is, tundra, mires, and forests, including hemi-, southern, middle, and northern boreal forests as well as sub-Arctic and alpine tundra. The areas were located mainly in protected nature reserves with low anthropogenic influence (Table S1 in Supporting Information S1).
Rásttigáisá study area is entirely above the treeline and has a mountainous and heterogeneous relief. Kilpisjärvi is mainly heterogenous mountain tundra but the lowest valleys dip into the mountain birch forests. Värriö is mainly boreal forest with open wetlands along gently sloping landscapes but the highest peaks ascend above the treeline. Tiilikka is dominated by peat bogs and the relief is relatively flat. Pisa is remarkably hilly and dominated by Norway spruce forests at lower elevations and Scots pine forests at hill tops. Hyytiälä is a mix of boreal forests and open wetlands and has relatively flat terrain. Karkali is a mix of broad-leaf forests (with temperate elements) and boreal needle-leaf forests in a relatively hilly terrain.

Measurement Sites Within the Study Areas
We conducted a random stratification to pre-select a suite of candidate measurement sites to maximally cover the main environmental gradients within the seven study areas (Aalto et al., 2022). This was conducted separately for each study area, except for Rásttigáisá study area and part of the Kilpisjärvi study area (see details below). The environmental variables that were used to stratify the environmental space vary from study area to another based on the area-specific characteristics but included variables such as total canopy cover, deciduous canopy cover, distance to forest edge, altitude, potential incoming solar radiation, and a topographic wetness index (SAGA wetness index) (Conrad et al., 2015).
In the random stratified site sampling methodology, first we masked the areas outside the study areas and other unsuitable areas such as lakes and extracted the remaining pixel information into a data frame. Next, we took a random subset of half of the remaining points (i.e., pixels) and used this data to reduce the multidimensional environmental space into its first three principal components with eSample function from iSDM R package Hattab & Lenoir, 2017). Then we took a sample of 100 points that maximally and systematically covered the shrunk environmental space. Because these 100 points may be located just next to each other, we repeated these procedures 100 times. Next, we pooled all the samples and used the frequency of selection for each point as the weight in the final random selection of points where we also kept a minimum distance of 100 m between the points. The final selection of the points, that is the measurement sites, was visually inspected from the histograms of the environmental variables. The final decision of each preselected measurement site was done in the field, for instance, by skipping sites that were logistically challenging or were too similar with other sites.
Selecting the measurement sites was carried out as described above, except for two study areas; Rásttigáisá study area and part of the Kilpisjärvi study area. All measurement sites in Rásttigáisá were based on a systematic study design , in which the sites were chosen using a stratified sampling based on elevation, potential incoming solar radiation, a topographic wetness index (SAGA wetness index) (Conrad et al., 2015), snow cover duration, and soil quality to represent main environmental gradients within the area. Some of the measurement sites (50/227 sites) in Kilpisjärvi were based on a systematic study design (Kemppinen et al., 2018;Tyystjärvi et al., 2021), which was complemented with sites in extreme soil moisture regimes (meltwater channels and ridge tops), snow conditions (short and long snow cover duration), and elevations (near mountain tops).

Soil Moisture Data
We measured soil moisture between 1 April 2020 and 30 September 2020 at all 503 measurement sites located within the seven study areas (Figure 1). We used soil moisture sensors (TMS-4 dataloggers; TOMST Ltd., Prague, Czech Republic), which measure soil moisture to a depth of c. 14 cm (Wild et al., 2019). We set the loggers to measure with a 15-min interval. The time-domain transmission method is used in the sensors to measure soil moisture (see Wild et al., 2019 for a detailed description of the sensor and the measurement technique). In short, the method is based on counting the number of electromagnetic pulses that travel within the counting unit within a unit time, and this number informs about the moisture content of the soil (high soil moisture decreases the number of pulses, low soil moisture increases them). The sensors produce raw time-domain transmission data, which we converted into volumetric water content (VWC%) using a conversion function adopted from Kopecký Note. Each measurement site was equipped with a soil moisture sensor. Climate data was derived for years 1981-2010 from the nearest meteorological stations (Table  S1 in Supporting Information S1).  (2021). We also tested conversion functions presented in Wild et al. (2019), which are specific to different soil types. For this purpose, we did a rough soil type classification in the field, where we assessed particle size of the mineral soil and measured the depth of the organic soil layer. However, we noticed that these soil type specific conversion functions resulted in highly unphysical VWC% values (e.g., <0% and >100%), especially in peatlands. Thus, we concluded that the conversion functions in Wild et al. (2019) are not applicable for the soil types present in our study areas, specifically, the organic soils in boreal peatlands. We tested the correlation of median July soil moisture across measurement sites between values calculated either with (a) the "universal" conversion functions or (b) the soil type specific conversion functions, and we found that the Spearman correlation coefficient was as high as 0.98. The difference between these two conversion approaches is in the absolute VWC% values. However, as the conversion did not greatly affect the relative order across the measurement sites, we decided to use the "universal" conversion function for all sites because it produces a much more plausible range of VWC% values.
Prior to the conversion, we removed all data where the raw soil moisture count was <200, as these counts are far outside the range and indicate that either the given sensor or its installation did not function properly (e.g., sensor was damaged or installed for instance into soils with air pockets). Next, we plotted all individual soil moisture time series month by month and inspected them visually. We identified all days when the data was clearly erroneous (e.g., times when the sensor was not in the soil or knocked down, for instance, by animals) and removed these days from the rest of the analyses. We also filtered out all soil moisture measurements of periods when the local soil temperature was below 1°C (measured with the same logger at the depth of 6.5 cm), because soil moisture values from frozen soils are invalid (Wild et al., 2019). In addition, we considered and tested additional quality checks, for example, to identify sudden but soon reversed drops in soil moisture, but after a careful inspection we concluded that the previous check-ups had already removed all suspicious data. For the analyses, we included only sensors with at least 90% temporal coverage (after the data cleaning explained above) during the measurement period calculated from the melting date (which varies from site to site; Figure 2) until the end of September.

Topographic Data
We used seven topography variables to explain the soil moisture patterns. These variables were calculated from a high-resolution digital elevation model (DEM). Geographic coordinates were recorded in the field using a high-accuracy Global Navigation Satellite System (GeoExplorer GeoXH 6000 Series; Trimble Inc., Sunnyvale, CA, USA) that provides centimeter-scale positioning accuracy under sufficient (clear-sky) conditions. The high accuracy of the measurement site locations enabled us to relate the soil moisture data with the DEM and its derivatives. The DEM was provided by the National Land Survey of Finland at 2-m resolution (see terrain maps in Figure 1) and was based on Light Detection and Ranging data covering entire Finland. A similar DEM product and its derivatives were not available for Rásttigáisá study area located in Norway, thus, we decided to exclude Rásttigáisá from the topographic analyses to ensure full comparability of the models.
We calculated a TWI (SAGA wetness index; hereafter SWI) (Conrad et al., 2015) using the SWI tool in the System for Automated Geoscientific Analyses software (SAGA GIS version 7.6.2). SWI is a multiple-flow-direction algorithm, which performs well as a proxy for soil moisture compared to other topographic wetness indices Riihimäki et al., 2021). We used the DEM at its finest resolution of 2 m, which performs well when calculating topographical wetness indices (Riihimäki et al., 2021, however, see also Ågren et al., 2014;Sorensen & Seibert, 2007). We used a filled DEM (following Wang & Liu, 2006), and 10 as the suction parameter (t) in calculating SWI.
We used a readily-available TWI product at 16-m resolution provided by the Natural Resources Institute Finland (Salmivaara, 2016). This TWI product is based on the algorithm by Beven and Kirkby (1979) and includes several improving pre-processing steps, such as the removal of the effect of roads blocking waterflow in the TWI calculations. The flow direction and flow accumulation rasters were calculated with the D∞ method (Tarboton, 1997).
We used a readily-available Cartographic depth-to-water index (DTW; Ågren et al., 2014;Murphy et al., 2009) product at 2-m resolution provided by the Natural Resources Institute of Finland (Salmivaara, 2020). This DTW product was based on the same DEM that we used in the SWI calculations. DTW was available with stream networks created based on various thresholds from 0.5 to 10 ha. After preliminary tests with the soil moisture data, we decided to use the DTW based on 0.5 ha threshold. , and the density plots are also presented in relation to each other in (d). The boxplots represent temporal variation in soil moisture within each area (e). In the box plots, the notches and hinges represent the 25th, 50th, and 75th percentiles, the whiskers represent the 95% percentile intervals, and the points represent the outliers. VWC%, volumetric water content.
We calculated topographic position indices (TPI) in SAGA-GIS software with two radii: 20 and 500 m (hereafter, TPI20 and TPI500). TPI is the relative elevational difference between the focal cell and the mean of its neighbors with a selected radius. Therefore, negative TPI values indicate a hollow and positive a hill or a ridge. Values close zero indicate even terrain or smooth slope.
We calculated the potential incoming solar radiation (hereafter, solar radiation) for the summer months (June, July, and August) to reflect the potential energy differences in surface energy balance between different locations. We used the same DEM as in previous calculations and conducted the calculations with the Potential Incoming Solar Radiation tool (Hengl & Reuter, 2008) in SAGA-GIS software. Solar radiation was first calculated monthly for the midmost days of each month with 1-hr intervals, and then, the three summer months were summed together. In the SAGA-GIS software, the algorithm takes into account the shadowing effect of the surrounding terrain as we calculated the visible sky (i.e., sky view factor) using a 10-km search radius and eight looking sectors and included this as an input in the solar radiation function alongside the DEM.
Lastly, we calculated a slope-penalized distance to the nearest waterbody index (hereafter, distance to water).
Here, we utilized all water features present in the Topographic database of Finland and the DEM provided by National Land Survey of Finland. We calculated the horizontal distance to the nearest water feature with the Accumulated Cost tool in SAGA-GIS software, and we set the slope as the cost surface in the calculations. We considered that the effect of the waterbody reaches further when the slope is even and decreases rapidly when slope is steep. We also considered that the effect of surface waters on soil moisture decreases rather rapidly as the distance increases between a water body and a measurement site. Consequently, we set all cost distance values >100-100. Then, we reversed the index so that index value 100 was given to (a) pixels which touch a water body and (b) pixels with zero slope and which are close to a water body. Conversely, index value 0 was given to pixels further away from water bodies over slopy terrain.

Soil Data
At each measurement point (n = 503), we measured organic layer depth and determined soil texture in the field. We determined the main soil type visually and by examining the soil between our fingers, and roughly classified the soils into five soil type categories, namely, clay, silt, sand, gravel, and organic soil (peat). We measured the depth using a thin metal rod (max. 80 cm). With the rod, we measured layers 0-10 cm to the nearest centimeter, and for layers >10 cm, we rounded the measurements to the nearest 5 cm.

Summary Statistics of Soil Moisture Variation
First, we characterized soil moisture and its variation in each study area by calculating mean soil moisture over all measurement sites and over the measurement period. We calculated average standard deviation (SD) both in time (over time within sensors in a given area) and space (over sensors in a given area). Furthermore, we explored how stable the spatial pattern of soil moisture was in the seven study areas by calculating Pearson correlation coefficient between all possible pairs of dates within the measurement period (Kachanoski & de Jong, 1988). All statistics were calculated on daily-aggregated (daily medians) soil moisture values.
Next, we explored how soil moisture mean was related to soil moisture variation following the methods in Brocca et al. (2012). Here, we calculated soil moisture mean over the measurement sites separately for each measurement date. Respectively, we calculated SD and coefficient of variation (CV) over measurement sites in each date. We did this for all study areas together and separately for each study area. Then, we inspected (a) how the mean and SD are related (hereafter, mean-SD relationship), and (b) how the mean and CV are related (mean-CV relationship), and here, we applied univariate linear regression models with either linear or quadratic terms. We tested the strength of the nonlinearity in the relationships by comparing the fits of the linear and polynomial models with the ANOVA. If the resulting p-value from ANOVA was sufficiently low (≤0.05) for the polynomial predictor, we deemed the relationship as nonlinear.

Statistical Models
We comprised three sets of models to analyze the influence of topography on spatial soil moisture patterns and its variability in time and across the study areas. We aggregated the soil moisture values into weekly averages for the 8 of 16 modeling. We considered only weeks when at least 66% of the study sites within a study area had full coverage of soil moisture data (i.e., excluding early spring when majority of the sites in a given area are still under snow).
First, we tested how the SWI, TWI, and DTW predicted the weekly mean soil moisture values in univariate linear models (LM) (Equations 1-3, Figures S1-S3 in Supporting Information S1). We chose these three predictors as the focus of the univariate analyses, as they are commonly used proxies for soil moisture (Ågren et al., 2014;Kopecký et al., 2021;Riihimäki et al., 2021). We fitted the LM in R (version 4.2.1; R Core Team, 2022) separately for each study area and week. We treated the weekly mean soil moisture as a response variable, and each of the three topography variables as an explanatory variable one at a time. We tested the predictive accuracy of the models with leave-one-out cross validation (LOOCV), in which each of the measurement sites were one by one removed from the data, rest of the sites used to fit the model. Finally, this model was used to predict soil moisture for the removed measurement sites. After all the measurement sites were predicted ones, these values were compared with the observed values by calculating a squared Spearman correlation coefficient (R 2 ).
Second, we fitted LM and generalized additive models (GAM) with six topography variables as predictors, namely, TWI, DTW, solar radiation, TPI20, TPI500, and distance to water, and one soil variable, namely, soil type as a categorical predictor (Equation 4; Table S2 in Supporting Information S1). The number of observations within soil types was highly unbalanced (e.g., very few measurement sites with clay as the main soil type). Thus, we aggregated the soil types into three categories, namely fine soil (including silt and clay), coarse soil (sand and gravel), and organic soil (peat). The soil type predictor represents the effect of soil structure on for instance water infiltration, but it may also the predicted soil moisture levels via controlling for the potential effect of the conversion from raw sensor data into VWC% values. Unlike LM, GAM allows non-linearity in responses. We fitted the GAMs with the restricted maximum likelihood estimation, and to avoid overfitting, we set the maximum dimension of the basis used to represent the smooth term as three. GAMs were fitted with the mgcv R package (Wood, 2011). Here, we tested how the influence of these predictors on weekly mean soil moisture patterns vary during the measurement period and across the study areas. SWI was excluded from these multiple regressions due to the high correlation between SWI and TWI. Compared to the SWI algorithm, the TWI algorithm is more commonly used, and the ready TWI product is openly available for entire Finland. The predictive accuracy of the models was tested with LOOCV R 2 and root mean squared error (RMSE). We calculated a permutation-based variable importance metric with the vi_permute function from vip package (Greenwell & Boehmke, 2020) to evaluate the relative importance of the predictors. To calculate the variable importance score, we compared the model fit (R 2 ) of a model fitted with the original dataset to a model in which, one at the time, each predictor was randomly permuted; if the permuted predictor is important, the R 2 will decrease greatly leading into a high importance value. The variable importance was calculated with 10 permutation rounds per predictor and study area. The method is model agnostic, and thus, can be compared across different modeling methods (Greenwell & Boehmke, 2020).
Soil moisture ∼ Topographic wetness index + Depth to water + Topographic position index (500 m radius) +Topographic position index (20 m radius) + Soil type + Distace to water + Solar radiation Third, we tested if the predictive performance of the multiple regression models was related to the overall wetness of the study areas. Here, we calculated the Pearson correlation coefficient between the weekly mean soil moisture and the R 2 value of the corresponding weekly models. We did this separately for each study area and modeling method. Next, we combined all these information into a single linear mixed effect model where we explained the R 2 value of the models by the weekly mean soil moisture. In the model, we included the modeling method (i.e., LM, GAM) as a categorical predictor and study area as random intercept (Equation 5). We fitted the model by using lme function from nlme R package (Pinheiro et al., 2022). We evaluated the significance of the relationships by using the F-test from anova function.
2 ∼ Mean soil moisture + Modelling method + (1|Study area) Finally, we tested how well the LM and GAM models fitted with the data from one study area predicted the soil moisture values of other study areas. Again, we used R 2 as a measure of predictive power. Here too, we used the mean weekly soil moisture as a predictive variable and the six topography variables as explanatory variables, but tested the model transferability only on weeks, in which the data from both the model fitting and prediction areas had at least 66% soil moisture data coverage. We tested the model transferability only across the study areas, and not in time (i.e., we did not test how model fitted for data from May performed when predicted to data from August, e.g.,).

Results
Soil moisture showed large spatial variation but often less pronounced temporal variations across the measurement sites within the seven study areas (Figures 1-3, Table 2). Measurement sites close to each other did not necessarily have similar mean soil moisture values (Figure 1). The soil moisture varied c. 0-60 VWC% at all study areas (Figures 1 and 2e). There were both slight drying and wetting trends in the study areas during the study period (Figures 2a and 2c). However, the spatial pattern remained relatively stable thorough the measurement period as the study areas showed high temporal stability (Table 2).
On average, there was more spatial variation within the study areas (10.9-16.7, SD across measurement sites) than temporal variation (2.9-6.4, SD within measurement sites) in soil moisture (Table 2). If the study areas were arranged by the average temporal variation within the areas, the study areas were nearly in their latitudinal order: northernmost site, Rásttigáisá (tundra), had on average the least amount of temporal variation, whereas the southernmost site, Karkali (hemiboreal), had the most ( Table 2). The amount of spatial variation in soil moisture within the study areas remained relatively stable during the measurement period.
The spatial soil moisture mean-SD soil relationships varied considerably across study areas (Figure 3a). We found little evidence for a unimodal mean-SD relation relationship since only one study area (Tiilikka) showed a significant unimodal relationship. However, when forcing the intercept of the models to zero (as in Brocca et al., 2012), the unimodal relationship was present. The mean-CV relationships followed the same non-linear relationship in all areas (except Tiilikka) where a decreasing trend levels out toward a minimum CV.
Modeling results showed that the soil moisture-topography relationships vary in time and across study areas. Weekly linear univariate models (Figure 4) showed that SWI had the highest overall predictive accuracy (averaged LOOCV R 2 over weekly models) in four of the study areas (Kilpisjärvi, mean R 2 = 0.26; Värriö, 0.41; Pisa, 0.39; Karkali, 0.38), TWI performed best in one area (Tiilikka, 0.43) and DTW in one area (Hyytiälä 0.26). DTW had typically the lowest predictive performance but the relative order of the three topographic varied and in some study areas (Pisa and Hyytiälä) in some individual weeks DTW scored the highest predictive performance. Also, the slope parameter estimates from the models varied greatly across study areas and also across weeks ( Figure S4 in Supporting Information S1).
Predictive performance of the multiple regression models showed similar pattern in space and time compared to the best performing univariate LM (GAM results in Figures 5a and 5b, LM results in Figures S5A and S5B in Supporting Information S1). Overall predictive performance was highest in models for Kilpisjärvi (mean LOOCV R 2 over weekly models and modeling methods = 0.42) followed by Värriö (0.39), Hyytiälä (0.33), Pisa (0.31), Tiilikka (0.30), and Karkali (0.18). The predictive performance was on average slightly higher for LM (0.31) than GAM (0.30), but the difference was not significant (paired Wilcoxon signed rank test, p = 0.53). On average, LOOCV RMSE was 11.8 for LM and 12.2 for GAM, and the difference was significant (p = 0.018). . Soil moisture mean-standard deviation (mean-SD) and mean-coefficient of variation (mean-CV) relationships. In (a), the colored lines represent the entirely data-driven quadratic mean-SD relationship, and the black lines represent relationships where the curve is forced to intersect the y-axis at 0. In (b), the colored lines represent the mean-CV relationship of the areas. The gray dots in the background represent the individual data points, that is, each dot represent one measurement date. Statistically significant non-linear relationships are marked as follows: *** p-value ≤ 0.001; ** p-value ≤ 0.01; * p-value ≤ 0.1.
The variable importance scores also varied across weeks for some areas (e.g., Karkali), while in others, the scores remained relatively stable (e.g., Pisa) (Figure 5c). On average, the most important variable for predicting soil moisture was TWI in three of the study areas (Tiilikka, Pisa, and Karkali), soil type in two areas (Kilpisjärvi and Hyytiälä) and the slope penalized distance to waterbodies in one area (Värriö). Variable importance scores for the topographic predictors followed similar patterns in GAM (Figure 5c) and LM ( Figure S5C in Supporting Information S1), but the scores varied greatly from area to another.
The predictive performance of the multiple regression model was weakly, positively related to the overall wetness of the soils in the study areas (Table  S3 in Supporting Information S1). Correlation coefficient was positive in 9 out of 12 occasions (modeling methods × study areas) but only three positive correlations were significant (p ≤ 0.05) and none of the negative relationships were significant. The GAM models for Karkali showed the highest correlation coefficient (0.74, p < 0.001), and the LM models for Kilpisjärvi the lowest (−0.31, p = 0.20). However, a mixed effect model which combined both modeling methods and all study areas, indicated a positive and significant (F-test, p > 0.001) relationship between the mean soil moisture and LOOCV R 2 .
Lastly, we tested how well the multiple regression models fitted in one study area performed when predicting soil moisture values of other areas ( Figures  S6 and S7 in Supporting Information S1). LMs had better model transferability on average (LOOCV R 2 = 0.23; LOOCV RMSE = 17.2) than GAMs (0.19; 19.6), and the difference was also significant (p < 0.001 in R 2 and p = 0.031 in RMSE). Models fitted with data from Kilpisjärvi (R 2 = 0.33; RMSE = 14.6) and Pisa (0.27; 14.0) performed best on average when these models were used for predicting soil moisture in other areas. Whereas models fitted with data from Hyytiälä (0.12; 20.2) and Tiilikka (0.16; 31.2) had low performance predicting soil moisture in other areas.

Soil Moisture in Boreal Forest and Tundra Environments
This study design and data recorded over several months revealed the locally heterogeneous patterns of soil moisture across and within boreal forest and tundra environments. Our results highlight that wide soil moisture gradients are present at small spatial extents, and that similar soil moisture gradients (0-60 VWC%) can be found locally from the tundra of northern Norway to the hemiboreal forests of southern Finland. Here, we also showed that soil moisture is relatively stable over time in the majority of the measurement sites (n = 503), and that across these environments, the spatial soil moisture pattern remains similar through the measurement period (April-September).
We found that there were also some measurement sites where soil moisture varied considerably in time. Studying these anomalous sites and their environmental conditions and ecosystem processes in more detail could be one important target for future studies to understand why and how they are different from the landscape matrices. Furthermore, if the annual precipitation increases but heat waves intensify in the northern environments due to climate change, the patterns of soil moisture variation can be very different in the upcoming decades (Samaniego et al., 2018). Thus, it is important to investigate if the heterogeneous soil moisture regimes can provide a buffer for ecosystems against increasing temperatures in the boreal forest and tundra  environments. Overall, there is a need for comprehensive study designs that monitor soil moisture across distinct environments and over large spatial and temporal extents (Dorigo et al., 2021).

General Patterns of Soil Moisture Variation
In general, the mean-SD relationship follows a unimodal pattern, in which sites with high and low mean soil moisture values express less variation compared to sites with intermediate mean soil moisture values (Pan & Peters-Lidard, 2008;Scaife et al., 2021). However here, we found only a weak signal of this unimodal mean-SD relationship. This is likely due to the temporal stability of soil moisture in the seven study areas. Overall, we found less temporal variation in mean soil moisture of the seven study areas in comparison to other studies (see e.g., Brocca et al., 2012;Dymond et al., 2021;Penna et al., 2009;Rosenbaum et al., 2012). Lawrence and Hornberger (2007) concluded that in humid regions, soil moisture variance decreases when mean soil moisture increases, and that in temperate regions, variance peaks at intermediate soil moisture contents. During the measurement period from April to September, we did not find extreme drying or wetting in the study areas. Therefore, Figure 5. Modeling results from generalized additive models (GAM) with six topography and one soil predictors. The lines represent the predictive performance of the weekly models as (a) the leave-one-out cross-validated (LOOCV) R 2 , and (b) LOOCV root-mean-square-error (RMSE). The thick lines represent the weekly results and the thin line the smoothed trends. In (c), the bars represent stacked variable importance of the predictors in the weekly GAM models. TPI, topographic position index; DTW, depth to water; TWI, topographic wetness index. Summary statistics of the predictors are provided in Table S2 in Supporting Information S1. the data do not cover situations where the very high or very low soil moisture values would force the mean-SD relationship into lower SD values in either ends of the mean soil moisture gradient. Thus, our findings are somewhat analogous with findings by Martínez-Fernández and Ceballos (2003). They found a monotonously increasing relationship between mean soil moisture and variance, but they found that mean soil moisture was always rather low, and thus, it likely covered only part of the whole potential soil moisture gradient.
We found a similar mean-CV relationship as in previous investigations (see e.g., Brocca et al., 2012;Famiglietti et al., 2008;Penna et al., 2009): a decreasing trend that eventually levels. Although, in comparison to previous investigations conducted in different ecosystems, we found much higher spatial variability across all measurement dates and thus across all mean soil moisture levels present in our data. Overall, we found that the spatial patterns remained relatively stable over time, and thus, these seven study areas are likely less prone to large temporal variation in SD across measurement sites. This is also in line with previous investigations, which have shown that in tundra environments, the spatial patterns of soil moisture are relatively stable during growing season (le Roux et al., 2013;Kemppinen et al., 2018;Tyystjärvi et al., 2021).

Soil Moisture-Topography Relationships
Soil moisture-topography relationships turned out to be context-dependent. This was expected, because of the other factors (e.g., vegetation, local temperatures) that control soil moisture and interact with topography. Thus, no single topography-based variable can manifest these relationships in all the study areas and conditions. We found that topography-based proxies cannot predict soil moisture patterns over the entire measurement period and across seven distinct study areas. This is important because topography-based proxies are widely used instead of field-quantified soil moisture in climate change and biodiversity modeling Riihimäki et al., 2021), although topography is only one of the drivers of soil moisture patterns (Albertson & Montaldo, 2003;Teuling, 2005;Wilson et al., 2004). Our results demonstrate why topography-based proxies for soil moisture should not be used without exhaustive investigations on their capability to characterize soil moisture patterns over space and time. This is especially important to consider in studies with large spatial extents and where high model transferability is required. By neglecting these investigations, topography-based proxies can lead to biased inferences on the importance of soil moisture, for instance, in ecological models (Kopecký & Čížková, 2010). Overall, we encourage careful consideration whenever soil moisture data are substituted with topography-based proxies Riihimäki et al., 2021) or optical proxies (von Oppen et al., 2021), regardless if soil moisture is modeled with statistical (Kemppinen et al., 2018) or process-based methods (Tyystjärvi et al., 2021).
Here, we showed that the influence of topography on soil moisture varies from study area and week to another. This means that, for instance, TWI is likely to perform well in heterogenous forest terrains (Tiilikka, Pisa) throughout the measurement period. Whereas in mountain tundra terrains (Kilpisjärvi), TWI tends to perform well only immediately after snowmelt when the landscape is very wet, and at other times, soil type was the most important predictor. This is in line with previous research which have reported temporal variation in the performance of topography as a proxy for field-quantified soil moisture (Ali et al., 2014;Riihimäki et al., 2021;Tague et al., 2010;Western et al., 1999), and here, we showed this with a harmonized study design across large geographic distances characterized by various environmental conditions, and for the entire snowless season of the boreal forest and tundra environments.
Our analyses did not show particularly clear trends in the strength of the soil moisture-topography relationships during the measurement period shared by all study areas, but there was a slight tendency that the topographic models had higher predictive accuracy early on the snowless season (i.e., after snowmelt) than later in the season.
In the sub-Arctic tundra, Riihimäki et al. (2021) also found that various TWI algorithms had stronger links to early season soil moisture (i.e., after snowmelt) than to the late season soil moisture. In other environments, the temporal variation of the soil moisture-TWI relationship has been linked to the overall wetness of the landscape (i.e., precipitation) (Ali et al., 2014;Western et al., 1999), and we too found that the model performance was slightly higher when the mean soil moisture across measurement sites was higher. Also, the temporal variation (or temporal instability) of the relationship is dependent on soil and vegetation factors (Coleman & Niemann, 2013).
Our results together with previous literature imply that in seasonally snow-covered environments, topography has stronger control on soil moisture shortly after snowmelt, and that the soil moisture-topography relationship weakens toward the end of the snowless season when other factors (e.g., evaporation, transpiration) play a greater role in controlling soil moisture patterns.
We found no single topography variable superior to others in predicting soil moisture patterns across all study areas, as the most influential predictors varied from area to another, and to some extent, also from week to week. For instance, in Karkali, the most relevant predictor changed from week after week, whereas in Pisa, one predictor (TWI) remained the most relevant for most of the measurement period. Furthermore, the transferability of the topographic models from one area to another varied. The models that were fitted with data from the study areas mostly covered by peatlands (Hyytiälä and Tiilikka) performed very poorly when they were used for predicting soil moisture in other areas. This indicates that soil moisture-topography relationships in wetland areas can be very different from other environments. Overall, our results highlight the need for detailed exploration and careful consideration of different topography variables in various environmental settings before applying them as proxies for soil moisture.

Challenges in Measuring Soil Moisture
Methodological challenges are various in soil moisture research, particularly related to field measurements (Robinson et al., 2008). Until recently, devices for continuous data have been expensive, which is impractical for detailed investigations over large spatial and temporal extents (Wild et al., 2019). Also, there are numerous device types that are based on different measurement techniques (Dobriyal et al., 2012;Romano, 2014;Yu et al., 2021). Moreover, different device types based on the same technique can provide significantly different results (Lekshmi et al., 2014;Walker et al., 2004), even the same device type can provide different values (Rosenbaum et al., 2010). In this study, we used only one device type to reduce uncertainty caused by different instruments (Wild et al., 2019), yet the manufacturer of the devices reports that error among the devices can be up to 5%. This is due to differences between the devices and soil types, and it cannot be controlled for, at least at the present moment. Soil type also influences the conversion of raw moisture measurements into VWC.
Here, we used a single conversion function adopted from Kopecký et al. (2021) as the previous functions (Wild et al., 2019) were created for different soils in the Czech Republic, and were considered difficult to apply for these northern study areas (typically high content of organic soil material) based on our preliminary tests (data not shown). As this specific sensor type is becoming more common among scientists across different ecosystems and soil types, we identify this as an important subject for development in the near future. Yet, at the same time, we consider that the data conversion is not a major issue in our study where the moisture gradient is wide and differences across measurement sites are considerable. Thus, with these approaches we managed to tackle critical challenges related to soil moisture measurements (Robinson et al., 2008).

Future Perspectives
One of main challenges remains for better understanding and modeling of spatio-temporal patterns of soil moisture and its effects on ecosystems: How to obtain field-quantified data on fine-scale patterns of soil moisture and its controlling factors with sufficient accuracy and extensive spatio-temporal coverage? Here, we found high temporal stability in soil moisture and its spatial patterns in boreal forests and the tundra of northern Europe, from peatlands to mountain tops. However, for instance in the tundra, soil moisture often varies greatly from one square meter to another (Kemppinen et al., 2019(Kemppinen et al., , 2021a. Consequently, further applications, such as ecological and biogeographical investigations, would require soil moisture data at a corresponding spatial resolution, and this is challenging regardless of methodology (e.g., statistical or process-based modeling). Here, we addressed the influence of fine-scale topography and soil properties on soil moisture, and the next step would be to incorporate the influence of fine-scale vegetation and local temperatures. Therefore, more data and models are needed on fine-scale soil moisture and its controlling factors, and particularly, more temporal data on vegetation and local temperatures. Our soil moisture data covered a single growing season and thus, it is likely that the data did not capture occasional extreme conditions (drought and flooding). Therefore, long-term measurements are needed to capture the full potential range in soil moisture dynamics. Finally, northern Europe is a seasonally snow-covered region, and thus, investigations on soil moisture-snow relationships are urgently needed for understanding the future soil moisture in boreal forest and tundra environments that are under rapid climate change.

Conclusions
Here, we present a rare case of intensive soil moisture investigations which cover large geographical and environmental gradients and the entire snowless season of the boreal forest and tundra environments in northern Europe.
We documented detailed soil moisture patterns at high temporal resolution in 503 locations within seven study areas from southern Finland to northern Norway. We found that soil moisture has high spatial variability in all seven study areas, and that the variability persisted the entire 6-month measurement period. We also found that the nature of soil moisture-topography relationships varied greatly across time and space. Overall, we highlight that wide soil moisture gradients can be present at small spatial extents, and that similar soil moisture gradients can be found locally across northern Europe from hemiboreal forests to the tundra.

Data Availability Statement
Data and code are openly available (Kemppinen et al., 2023). Field Station for their support during fieldwork. We acknowledge funding for fieldwork and equipment by the Nordenskiöld samfundet, Tiina and Antti Herlin foundation, and Maa-ja vesitekniikan tuki ry. JK was funded by the Academy of Finland (project number 349606), and the Arctic Interactions at the University of Oulu and Academy of Finland (project number 318930, Profi 4). PN was funded by the Academy of Finland (project number 347558), the Nessling foundation, and the Kone Foundation. TR was funded by the Doctoral Program in Geosciences, University of Helsinki. JA acknowledges the Academy of Finland Flagship funding (project number 337552). Permission to carry out fieldwork was granted by Metsähallitus.