Ongoing Fen–Bog Transition in a Boreal Aapa Mire Inferred from Repeated Field Sampling, Aerial Images, and Landsat Data

Northern aapa mire complexes are characterized by patterned fens with flarks (wet fen surfaces) and bog zone margins with Sphagnum moss cover. Evidence exists of a recent increase in Sphagnum over fens that can alter ecosystem functions. Contrast between flarks and Sphagnum moss cover may enable remote sensing of these changes with satellite proxies. We explored recent changes in hydro-morphological patterns and vegetation in a south-boreal aapa mire in Finland and tested the performance of Landsat bands and indices in detecting Sphagnum increase in aapa mires. We combined aerial image analysis and vegetation survey, repeated after 60 years, to support Landsat satellite image analysis. Aerial image analysis revealed a decrease in flark area by 46% between 1947 and 2019. Repeated survey showed increase in Sphagnum mosses (S. pulchrum, S. papillosum) and deep-rooted vascular plants (Menyanthes trifoliata, Carex rostrata). A supervised classification of high-resolution UAV image recognized the legacy of infilled flarks in the patterning of Sphagnum carpets. Among Landsat variables, all separate spectral bands, the Green Difference Vegetation Index (GDVI), and the Automated Water Extraction Index (AWEI) correlated with the flark area. Between 1985 and 2020, near-infrared (NIR) and GDVI increased in the central flark area, and AWEI decreased throughout the mire area. In aapa mire complexes, flark fen and Sphagnum bog zones have contrasting Landsat NIR reflectance, and NIR band is suggested for monitoring changes in flarks. The observed increase in Sphagnum mosses supports the interpretation of ongoing fen–bog transitions in Northern European aapa mires, indicating significant ecosystem-scale changes.


INTRODUCTION
Climate-driven changes are likely in northern mires (peat-accumulating wetlands), as many mire types depend on specific climatic conditions (Ruuhijä rvi 1960;Rydin and others 1999;Kuznetsov 2003;Luoto and others 2004). Climate warming may cause a fall in the mire water level (Gong and others 2012), and changes in plant communities (Kokkonen and others 2019;Kolari and others 2021) and carbon cycling (Mä kiranta and others 2018). Some studies have reported drying in mires (Swindles and others 2019; van Bellen and others 2018) and projected a decreasing carbon sink in future (Leifeld and others 2019; Chaudhary and others 2020; Hopple and others 2020; Hugelius and others 2020), while other studies have indicated increase in carbon sink due to increasing productivity (Charman and others 2013; Gallego-Sala and others 2018). Simultaneously, mires are threatened by peat extraction, clearing for cultivation, and forestry drainage, making them the most threatened ecosystems in Europe (Janssen and others 2016).
Although climate and land use impacts on mires have gained interest, studies have only recently focused on changes in boreal mire ecosystems, such as aapa mires (Sallinen and others 2019; Granlund and others 2021; Kolari and others 2021; Robitaille and others 2021; Primeau and Garneau 2021). Aapa mires are characterized by wet central fens and Sphagnum-dominated bog margins (Ruuhijä rvi 1960;others 2007, 2017;Sallinen and others 2019). The fens have mosaic patterns of sparsely vegetated wet depressions called ''flarks'' and drier hummock strings. These features are common to patterned fens covering wide boreal regions in North America (Glaser and others 1981;Foster and King 1984) and Russia (Yurkovskaya 2012). Some evidence exists of a decrease in flark area and an increase in bog vegetation in patterned fens in recent decades (Tahvanainen 2011; Granlund and others 2021; Kolari and others 2021; Robitaille and others 2021). The transition from a minerotrophic fen to a Sphagnum bog indicates an ecosystem-level change and, most importantly, a potential increase in carbon accumulation (Tahvanainen 2011;Loisel and Yu 2013;Robitaille and others 2021) and a decrease in methane emission (Zhang and others 2021).
Because aapa mires cover vast areas, remote sensing methods are needed for effective monitoring. Aerial image analyses provide efficient means to estimate land cover in mire ecosystems (Ihse 2007;Waser and others 2008;Tahvanainen 2011;Edvardsson and others 2015). Gray-tones of historical aerial images connect to hydro-morphological patterns of aapa mires; wet flarks appear blackish, while Sphagnum lawns have whitish tones (Tahvanainen 2011). Furthermore, visual interpretations of aapa mire patterns are relatively easy in the case of individual aerial images. However, processing a long time series of high comparability is challenging, since the availability, accessibility, and quality of old aerial images vary considerably.
Satellite images offer certain benefits over aerial images. In particular, Landsat images have global coverage, are freely available, and cover a multidecadal timespan that allows versatile applications for monitoring ecosystems (Zhu 2017), including boreal wetlands (Terentieva and others 2016;Wulder and others . 2018;Sirin and others 2020). Among satellite-derived spectral variables, the Normalized Difference Vegetation Index (NDVI; Rouse and others 1973) is among the most frequently used globally and in northern mires (Torabi Haghighi and others 2018; McPartland and others 2019). NDVI reflects vegetation abundance and greenness, but it has also been connected to soil moisture (Atkinson and Treitz 2013) and drying of degraded bogs (Š imanauskien _ e and others 2019). Many other indices have also been developed to identify vegetation characteristics and other landscape properties (Li and others 2013;Xue and Su 2017;Ma and others 2019).
Field and laboratory studies have attempted to distinguish Sphagnum mosses from vascular plants based on spectral features (Bubier and others 1997;Pang and others 2020) and to connect spectral signals of Sphagnum to moisture content, desiccation damage, and photosynthetic activity (Harris 2008;Meingast and others 2014;Lees and others 2019). Moreover, spectral features of Sphagnum differ from vascular plants in visible, near-infrared (NIR), and short-wave infrared (SWIR) regions (Whiting 1994;Bubier and others 1997;Harris and others 2005; Lees and others 2020; Pang and others 2020), while water absorbs NIR radiance effectively (Jones and Vaughan 2010;Lillesand and others Ongoing Fen-Bog Transition 2015; Mondejar and Tongco 2019). Spectral indices derived from NIR water absorption were suggested for hydrological monitoring of raised bogs (Harris and others 2006;Harris and Bryant 2009), based on a correlation with Sphagnum moisture content (Harris and others 2005). In another case, a SWIRbased model was suggested for water-table monitoring (Burdun and others 2020). Middleton and others (2012) used airborne hyperspectral data and demonstrated that the NIR region provided the most information for mire type separation. To explain variation in aapa mire vegetation, NIR, red, and green bands and NDVI performed well among combined spectral and topographic data (Rä sä nen and others 2019, 2020), indicating potential utility in the detection of vegetation changes in aapa mires.
In this study, we combined analysis of Landsat satellite data, aerial image time series, and repeated field sampling to explore recent vegetation changes in a south-boreal aapa mire complex. The study site was ideal for testing remote sensing methods in detecting Sphagnum increases in aapa mires, as it has a characteristic pattern of central flark fen with bog zones at the margins, and accurate field data of vegetation and hydro-topographic patterns from 1959 (Tolonen 1959). We repeated the field study in 2018, supported by ultra-high-resolution UAV (Unmanned Aerial Vehicle) imaging, to provide ground validation for coarser-scale remote sensing. We test the correlation between fen area changes derived from aerial images (1947-2019) and from Landsat data. We hypothesized that recent warming and catchment disturbance have accelerated succession toward bog vegetation and that the change is reflected in remote sensing data. The study area's location near the southern distribution limit of the aapa mires suggests particular sensitivity to climate warming.
The annual average temperature is + 3.5°C, and the annual precipitation is 552 mm (1961, ClimGrid data, Finnish Meteorological Institute, described in Aalto and others 2016. From 1961 to 2018, annual temperatures rose ca. 0.3°C and annual precipitation ca. 13 mm per decade. The studied mire is largely undrained, but it is surrounded by ditches ( Figure 2). According to old aerial images, ditches in the eastern margin existed in 1947. Since then, more drainage has taken place, and at present, all mire margins are ditched.  Figure 1. Location of the study area and the distribution of mire complex types in Finland. The points represent peatland areas with undrained area of at least 50 hectares according to Finnish Environment Institute's GIS and aerial photograph survey (see Sallinen and others 2019). Climatic mire zones according to Ruuhijä rvi (1960Ruuhijä rvi ( , 1988. Base map: Helcom.

Remapping of the Flark Fen Special Study Area
The study site was first studied in detail in 1959 by Professor Kimmo Tolonen (Tolonen 1959, see Supplement 1.1). The data included a flark fen special area (0.5 ha) in the wettest minerotrophic part of the mire, with graphical descriptions of vegetation with remarks on dominant plant species and delineations of open water flark patterns. In addition, vegetation cover was estimated from 80 vegetation plots (20 9 20 cm) along a 16-m transect marked in the same graph. We visited the study area in July 2018 and repeated the vegetation sampling. We collected water samples (n = 3) that were analyzed for pH and several mineral elements by inductively coupled plasma-mass spectrometry (ICP-MS), of which we report calcium concentrations. Prior to field work, we located the flark fen special area in an aerial image from 1963 and georeferenced the transect and the flark patterns based on 1959 field work, with some corrections to the original drawing. In the field, we located the 1959 flark boundaries and the transect with a GNSS RTK device (TOPCON).
The starting point of the transect was on an S. fuscum hummock, which enabled the location of the point in the field. Since the precise starting point was still somewhat uncertain, we collected new data from two contiguous transects with 1-m intervals (transects A0 and A1) to consider the possible inaccuracy of the location. Along both transects, the percent coverage of each plant species, open water, and mud-bottom surface were visually estimated, and these percentages were later converted to the scale used by Tolonen in 1959 (Supplement 1.1).
The main directions of variation in the plant community composition along the transects were analyzed with detrended correspondence analysis (DCA; Hill and Gauch 1980), using combined 1959 and 2018 vegetation plot data. DCA is an indirect ordination method, that is, it is used to reveal vegetation compositional gradients, which can later be interpreted in relation to environmental variation. Data from transects A0 and A1 were separately combined as paired observations with the 1959 data and then subjected to DCA. Cover values for ''open water'' and ''mud-bottom'' were included because no plants occurred in four plots in 1959. The species S. majus and S. jensenii were treated jointly as they were treated as S. dusenii in the 1959 data, and Cladonia lichens were treated at the genus level. In addition, S. parvifolium (syn. S. angustifolium) of Tolonen (1959) was interpreted as S. balticum, as correct species identification in 1959 could not be ascertained. In 2018, S. balticum was present in the plots where S. parvifolium was recorded in 1959. Similarly, S. recurvum (unknown in Finland, unclear taxon) of Tolonen (1959) was interpreted as S. pulchrum. These adjustments made the analyses conservative against false interpretations of changes due to misidentification or nomenclature. DCA was performed in PC-ORD (version 7.04, McCune and Mefford 2016).

UAV Image Interpretation
We obtained an ultra-high-resolution UAV image in July 2018 to illustrate present-day vegetation patterning in the special study area and to aid interpretation of aerial and satellite images. The drone used was an eBee Plus RTK, and the RGB image had a spatial resolution of 3.6 cm. After creating an orthomosaic, supervised image classification was performed in ArcMap 10.5. The goal of the classification was to visualize modern mire surface patterning and to determine whether the flark patterns of 1959 could still be recognized in the landscape. The selection of training samples for image classification was based on field notes from July 2018 and visual interpretation of the image. Five final classes were identified: 1) open water, 2) Carex-Menyanthes communities, 3) Sphagnum carpets, 4) S. papillosum lawns, and 5) Sphagnum hummocks. Sphagnum carpets consisted mainly of flark-level Sphagnum spp, (S. pulchrum, S. jensenii, and S. lindbergii). Lawns were characterized by S. papillosum, but S. compactum and Carex rostrata were also abundant, along with abundant sedge litter, which was difficult to distinguish from S. papillosum in the drone image. Hummocks were mainly composed of S. fuscum and S. rubellum. Each class was represented by 10 samples of more than 50 pixels which we found sufficient given the small size of the area of interest (0.5 ha). The selected classifier was the support vector machine (SVM) and a majority filter (3 9 3 window) was applied to remove the salt-and-pepper effect of the classification results.

Aerial Image Analysis
We explored a time series of aerial images to study changes in the flark area throughout the entire mire. The focus was on the flark fen special study area of Tolonen (1959), of which we had detailed field data from 1959 and 2018. We used ArcGIS 10.7.1 for this and the subsequent remote sensing and GIS operations, if not mentioned otherwise. First, we delineated the mire catchment and the aapa ecosystem under consideration, based on a digital elevation model (DEM KM2) of the National Land Survey of Finland (NLSF), originally with 2-m horizontal resolution but averaged to a 10-m resolution. We used the GIS method described by Sallinen and others (2019), with the idea of automatically defining the topographic catchment for an outlet area in the mire from where surface water flow exits. Open mire was extracted from the forested and treed areas in the catchment using a peatland data set (NLSF open data) and vegetation height data (resolution 2 m, refined in the Finnish Environment Institute SYKE from LIDAR data of NLSF).
Archival black-and-white aerial images (dates 1947-07-17, 1988-05-20, 1997-06-07) were purchased from NLSF with differing resolutions (0.33-0.48 m). We resampled the images to 0.5-m resolution, as in the newest digital aerial image used (2019-06-06). A low-pass filter (mean filter with a 3 9 3 window) was used to remove the noise present in the old images. We applied supervised classification (maximum likelihood) to the images representing only open mire areas. Based on visual assessment, we selected training samples for two classes, flark and non-flark, with approximately one hundred pixels for both in each image. We performed the classification multiple times with different training samples and selected the results that best expressed the visually observable flark boundaries. Tolonen's (1959) open water flark delineation of the special study area, field notes, and the drone image from 2018 provided ground validation data. Particular attention was given to the consistency of flark boundaries in different images. The flark classification results were subsequently converted to polygon format, and misclassified features, like tree shadows or S. fuscum hummocks incorrectly classified as flarks, were corrected manually. From the resulting layers, we calculated flark areas and their temporal changes.

Satellite Data Analysis
We used the USGS Earth Explorer to search cloudless Landsat images acquired as close to the dates of the aerial images as possible. The dates of the selected Landsat images were 1988-05-12, 1997-06-13, and 2019-06-19. In addition, we selected 18 Landsat images acquired between 1985 and 2020. To minimize the phenological differences between images, these were acquired in July or August (Rafstedt and Andersson 1982;Saarinen and others 2018). The product identifiers of all Landsat images used in this study are presented in Supplement 2.1. Among suitable satellites, Sentinel-2 or SPOT would have provided higher spatial resolution data, but Landsat covers a longer period for the assessment of decadal-scale ecosystem changes.
We used the Collection 1 Tier 1 products, which are georegistered within £ 12 m radial rootmean-square error (USGS 2016). Nevertheless, the location accuracy of all images was checked. For preprocessing, we used the R function calcAtmosCorr in the ''Satellite'' package (Nauss and others 2019). The function computes an atmospheric scattering correction and converts the sensors' digital numbers to reflectance values using absolute radiance correction and the dark object subtraction method DOS2 (Chavez 1996). Our data comprised images from the three satellites Landsat 5 (TM), 7 (ETM +), and 8 (OLI), exhibiting differences in their technicalities. Most importantly, the near-infrared (NIR) band is narrower in Landsat 8 than in TM and ETM + (USGS 2020). To improve the continuity between the images, we applied the corrections proposed by Roy and others (2016) ( Table 1).
We created a 30-m cell size grid layer to match Landsat pixels (30 9 30 m), and based on the aerial image classifications described above, we calculated the flark area of each pixel in different time phases (1988, 1997, and 2019). Subsequently, we retrieved surface reflectance data from the Landsat images of the same time phase and calculated several vegetation and water indices from the data ( Table 2). We tested the correlation between flark area per pixel and the spectral bands and indices with the non-parametric Spearman's rank correlation test because the assumption of a normal distribution was neither met with the flark area nor with most of the Landsat-derived variables, as judged by the visual assessment of histograms and the Shapiro-Wilk test of normality. The variables that showed the strongest correlations consistently in all time periods were selected for a more detailed examination of spatial and temporal variation, using the time series of Landsat images from 1985 to 2020. We also examined variation in NDVI, as it is the most used spectral index and was previously suggested for peatland monitoring (for example, Š imanauskien _ e and others 2019). We compared the selected spectral variables between the flark fen special area and the surrounding open Sphagnum bog in the studied mire, as well as reference sites, a nearby raised bog, and a coniferous forest area in a nature reserve located in the same Landsat images.

Weather Data
We obtained data on daily temperature and precipitation from the Finnish Meteorological Institute (open data) and calculated mean temperature and precipitation sum for 2, 5, 10, and 30 days before the acquisition dates of Landsat images used in the satellite analyses. For 1985-2018 we used a gridded, permutation based ClimGrid data set (Aalto and others 2016), and for 2019 and 2020, data from the closest meteorological weather station (Pelmaa, Seinä joki). We determined correlations between weather variables and residual variation from linear trends in the time series of Landsat variables. The purpose of this analysis was to test if shortterm weather variation would explain residual variation of long-term trends of Landsat data. Weather variation is autocorrelated with the longterm time series due to climate warming and with seasonal variation due to varying day of year of Landsat data. The data are not sufficient for more detailed account of weather impacts on Landsat inference.

Vegetation Changes in the Flark Fen Special Study Area
In 1959, the flark fen special area was characterized by large open water flarks alternating with Sphagnum papillosum lawns ( Figure 3, Supplement 1.1). Dominant vegetation types were S. papillosum carpets, short-sedge S. papillosum lawn fens, and S. papillosum-Carex rostrata communities. S. cuspidatum was found in one flark, while others had open water and flark-level vascular plant species, namely Menyanthes trifoliata, Carex rostrata, and C. limosa that were growing mainly on the flark margins. Warnstorfia fluitans, S. pulchrum, S. lindbergii, Carex canescens, Trichophorum cespitosum, and Rhynconphora alba occurred sparsely in the special study area.
Since 1959, Sphagnum mosses increased in the flarks, and in 2018, only small remnants of the formerly large open water pools existed ( Figure 3). Among flark-level species, M. trifoliata and C. rostrata were still abundant in 2018, and in some parts, they had even increased, whereas W. fluitans and S. cuspidatum were not observed. The most noticeable change was the increase in S. pulchrum that expanded over the flarks, in addition to S. papillosum and S. compactum. Two of the Sphagnum  (Figure 3).
The decline in open water area, expansion of Sphagnum carpets (mainly S. pulchrum, S. jensenii, and S. lindbergii), and the increase in flark-level vascular plants (particularly C. rostrata and Menyanthes trifoliata) were also observed with supervised classification of a high-resolution drone image (Figure 4). In the classified image, the connected pattern of the former open water flarks reappeared mainly as a Sphagnum carpet, differing from surrounding S. papillosum lawns.
The water samples collected in the flark fen special study area had pH values of 3.71, 4.02, and 4.08 and calcium concentrations of 1.93, 2.29, and 3.97 mg/l. Concentrations of other main elements are reported in Supplement 1.2.

Vegetation Plot Ordination
In 2018, vegetation plot data was collected from two transects (A0 and A1) that were separately combined with the 1959 plot data and subjected to DCA. Both analyses resulted in similar plot positions in the ordination space ( Figure 5). In the twodimensional DCA ordinations, the first axis was interpreted as a gradient in water-table depth. Common bog hummock species, such as S. fuscum, Calluna vulgaris, and Empetrum nigrum, characterized the plots on the left of the ordination. Plots characterized by flark-level species Warnstorfia fluitans, Menyanthes trifoliata, Carex limosa, and Carex rostrata, as well as ''open water'' and ''mud-bottom,'' spread to the right in the ordination, and were only represented by the 1959 data.

Aerial Image and Landsat Analysis
According to aerial image analysis, the flark area decreased over the whole mire area by 46% from 1947 to 2019. In the special study area, the flark area decreased by 76% (Figure 6). The aerial image-derived flark area had the strongest Landsatvariable correlations (r S between -0.84 and -0.89) with green, red, and NIR bands, and among the indices, with GDVI in the 1988 and 1997 data (Table 3). In the 2019 data, green (r S = -0.80) and red (r S = -0.81) bands' correlations remained high, while correlations of NIR (r S = -0.53) and GDVI (r S = -0.31) were weaker but still consistent and significant. Short-wave infrared (SWIR1, SWIR2) and blue bands also had consistent negative correlations with flark area (r S between -0.51 and -0.79). The automated water extraction index (AWEI) had a weak but consistent and significant positive correlation with flark area (r S = 0.67, 0.39, and 0.48 for 1988, 1997, and 2019, respectively). Among the rest of the Landsat-indices, the difference vegetation index (DVI) and green vegetation index (GVI) had weaker correlations, but they were still consistent between the studied years. All other indices either lacked significant correlation with the flark area or had disparate correlations between years (Table 3) still visible in the red band and AWEI but not in the NIR band and GDVI. In the case of NDVI, the difference between the central flark fen and the surrounding areas became apparent in 2020, with higher index values in the flark fen. Between 1985 and 2020, NIR and GDVI increased in the central flark fen and decreased in the surrounding bog area. The spatiality of NDVI changes was largely similar, with the difference that the NDVI also increased in the northeastern corner. AWEI declined throughout almost the entire study area. Red re-flectance increased in the western part of the study area and decreased in the north and northeast, and the changes did not appear to be related to the changes in the flarks. This finding also concerns blue, green, SWIR1, and SWIR2 reflectance (Supplement 2.2). The series of 18 Landsat images revealed changes in mean NIR and red bands, and NDVI values in different parts of the Mahlaneva mire and nearby reference sites for a 35-year period (Figure 8). The open Sphagnum mire and reference raised bog had the highest NIR reflectance among the studied areas. From the beginning of the time series until 2006, the flark fen had a distinctly lower NIR reflectance, similar to reference coniferous forest, but from 2006 to 2020, the NIR reflectance of the flark fen rose near to the level of the reference raised bog. The open Sphagnum mire and both reference sites had no significant trends of NIR reflectance. The NDVI showed synchronous fluctuations between sites (low values 2006-2010) but no significant trends. The NDVI level of the flark fen did not differ from that of the Sphagnum bogs. Only the reference coniferous forest differed from the other sites with higher NDVI values throughout the study period.
Weather data did not show consistent correlations with residual variations in the linear regression fit to temporal trend of NIR or red bands or NDVI (p > 0.1, n = 18), but day of the year had a significant positive correlation (r = 0.667, p = 0.002). Together, the year and day of the year of Landsat images explained over 80% of the variation in NIR reflectance in the flark fen when using linear regression (r 2 = 0.818, p < 0.001). For

Axis 1
Axis 2  instance, the low NIR levels in 1998, 2010, and 2020 had early timing compared to preceding and following time points (Figure 8).

Recent Changes in Vegetation and Hydromorphology
The decrease in flark area was detected in the studied aapa mire complex by repeated field sam- trifoliata and Carex rostrata, also increased. The spatial patterning of present-day vegetation closely conformed to the shapes of the former open water flarks. These results gave reliable ground validation to the Landsat analysis, according to which the NIR surface reflectance of flark fen was originally low but rose to the same level as the surrounding bog and the reference raised bog site, reflecting the shift toward Sphagnum-dominated vegetation.
Our results indicated two simultaneous processes: a diffuse increase in Sphagnum cover over most of the flark fen area and a narrowing of the continuous central flark fen zone ( Figure 6). The latter process may be interpreted as lateral expansion of bog zone margins that widely exceeded 20 m during the study period in the most clearly altered places along the western side of the flark fen. This translates to a lateral expansion rate of ca. 30 cm y -1 .
Natural succession from minerotrophic flark fen toward Sphagnum bog vegetation has probably been ongoing already before the first sampling in 1959, as concluded by Tolonen (1959) and more widely suggested by Tolonen (1967) in eastern Finland. However, our results from the NIR time series indicated that an increasingly progressive phase of this development occurred only after 2006. In contrast, aerial image analysis showed a marked decrease in the flark area between 1988 and 1997. The decline in water index AWEI throughout the mire between 1985 and 2020 must also indicate overgrowth of the water surface by Sphagnum carpets but not lowering of water level, since infilling of flarks does not equate to drying, as the moss carpets are essentially floating rafts forming ''quaking'' mire surface.
Leveling performed in Tolonen's (1959) study did not reach the mineral soil, which could have confirmed that the water level had not descended. However, we found indirect evidence. According to the author of the 1959 data, K. Tolonen, the flark pool, across which the vegetation plot transect was  Ongoing Fen-Bog Transition located, was not nearly as deep in 1959 as it was during 2018 field work (see photographs in Supplement 1.4). Second, the height difference between the leveling points of the bog margin and central fen (cross-transect X, Supplement 1.5) had leveled out from 46 cm in 1959 to 16 cm in 2018, suggesting that the peat surface and water level of the central fen area had ascended by 30 cm. However, the bog margin surface elevation may have also changed, for example, descended by subsidence or raised by peat accumulation. Several species that had increased in the flark fen of our study site (S. pulchrum, S. papillosum, Menyanthes trifoliata, Carex rostrata) are indicators of weak minerotrophy (Eurola and others 2015). The endangered species S. molle, which was recorded in the 2018 field survey, is also a minerotrophic indicator species (Heikkilä and Lindholm 1988). We interpret that the flark fen was still weakly minerotrophic, although some indicators were deep-rooted species that may survive prolonged periods after ombrotrophication (Gorham 1950(Gorham , 1957. The water samples also indicated minerotrophy, showing calcium concentrations between 1.9 -3.9 mg/l in the special study area, similar to fens elsewhere in Finland (Tahvanainen and others. 2002;Tahvanainen 2004). Interestingly, the mud-bottom hollow marked in the 1959 graph had been vegetated by Sphagnum mosses and Rhynchospora alba (Supplement 1.1). We also found another mud-bottom hollow with bare peat surface, located within the 1959 graph area but not noted by Tolonen (1959) (Figure 3). The formation of mud-bottom hollows is known to be highly dynamic, and they are often inundated during wet seasons and become rapidly overgrown by Sphag-num carpets, although larger mud-bottom hollows may persist over decades (Karofeld and others 2015).

Interpretation of the Causes of the Observed Changes
The observed infilling of open water flarks by Sphagnum carpets and the lateral expansion of Sphagnum-dominated bog zones follow the general bog succession scheme (Gorham 1957;Kuhry and others 1993;Belyea and Clymo 2001). Autogenic development from a minerotrophic fen to an ombrotrophic Sphagnum-dominated bog is caused by long-term peat accumulation and gradual isolation from minerogenic water level, and consequent vegetation changes (Gorham 1957;Kuhry and others 1993;Belyea and Clymo 2001). In our study area, such a development process had already started before drainage in the catchment or the first field survey in 1959. In his thesis, Tolonen (1959) interpreted the flark fen special study area as ''the last battleground of aapa mire and raised bog'' (Supplement 1.1). His cross-sectional peat profiles revealed a steady narrowing of Carex peat strata until recent times, indicating that minerotrophic fen vegetation had likely been decreasing since the sub-boreal period (Supplement 1.6). The narrowing of the fen by the lateral expansion of Sphagnum peat from the margins toward the center continued until recently. Our survey in 2018 demonstrated Tolonen's (1959) prediction of declining fen habitat. The gross topography of our study area is still slightly concave but already quite close to horizontal (Supplement 1.5), thus approaching the transitional phase that would be expected in the case of autogenic succession resulting from peat accumulation.
Although the development from fen to bog vegetation ultimately results from autogenic succession, this process may have been accelerated by allogenic factors (Hughes 2000;Hughes and Barber 2004), such as catchment drainage (Tahvanainen 2011) and recent warming (Kolari and others 2021). In our study area, the mean annual temperature has risen ca. 0.3°C, and the annual precipitation sum has increased ca. 13 mm per decade since the early 1960s (Finnish Meteorological Institute, open data), which is similar to general climate trend in Finland (Aalto and others 2016) and continues the longer period of warming after the Little Ice Age (LIA) (Luoto and Nevalainen 2017). The effective temperature sum has risen on average from 1100°C (1961-1979)   that has been considered to divide the climatic ecotone between raised bog zone to the south and aapa mire zone to the north in Finland (Ruuhijä rvi 1960). Sphagnum productivity has been shown to increase with similar warming (Dorrepaal and others 2004;Gunnarsson 2005;Breeuwer and others 2008), and recent warming may have accelerated the expansion of Sphagnum carpets in the flarks. Sphagnum mosses also benefit from decreasing snow cover (Zhao and others 2017), as well as from longer growing seasons, if the moisture supply is not limited (Loisel and others 2012;Kü ttim and others 2020;Bengtsson and others 2021). Moisture conditions in our study area have remained well sufficient for Sphagnum mosses to take advantage of warming. A recent increase in Sphagnum mosses over fen vegetation, coincident with a 20-year period of significant warming, was also observed in a pristine rich fen in eastern Finland without any drainage in the catchment (Kolari and others 2021). Paleoecological reconstructions have also indicated fen-bog transitions in patterned fens after LIA or coinciding with twentieth century warming in North America (Loisel and Yu 2013;Magnan and others 2018;van Bellen and others 2018;Primeau and Garneau 2021;Robitaille and others 2021), while our study adds one case of this phenomenon in Northern Europe.
In addition to the direct effects of warming, it is possible that water level fluctuations have had impact on vegetation (Granath and others 2010). A detailed account on possible water level effects is not possible, however, as no monitoring existed. Drainage activities in the mire catchment may have reduced the input of minerogenous water and affected water quality and flow rate even when the water level was not affected. Similar situation of margin drainage was documented in the Valkeasuo aapa mire in eastern Finland, where catchment disturbance caused hydrology-mediated changes in mire vegetation in an allogenic process (Tahvanainen 2011). Like in our case, the central mire area remained wet and deep-rooted minerotrophic vascular plants were present, while open water flarks were invaded by quaking Sphagnum moss carpets. In our interpretation, recent warming and catchment drainage have possibly contributed to driving the vegetation changes in our study site. As aapa mire pools and open water flarks are of secondary origin (Foster and others 1988;Kutenkov and Philippov 2019), their observed infilling by Sphagnum carpets can be considered secondary hydroseral development.

Performance of Remote Sensing in the Detection of Aapa Mire Changes
The comparison of old and new aerial images is an intelligible way of detecting changes in mires. Previously, historical aerial images have been successfully used in studies examining the thawing of permafrost in palsa mires, for example, as the spatial resolution enables the detection of small-scale changes (Pironkova 2017;Zuidhoff and Kolstrup 2000). In aapa mires, the contrast between wet flarks and Sphagnum covered areas makes changes readily visible, but experience is required on interpretation and evaluation may be misguided, for example, by the temporal fluctuation of water level (Ihse 2007;Tahvanainen 2011). Measuring areal changes is yet more difficult, since it requires the localization of flark boundaries that can be diffuse in nature, and varying image quality adds to the difficulty of classifying old images (Thibault and Payette 2009;Pironkova 2017). For meaningful comparisons, the determination of flarks must be repeatable between different images. This can be verified with good ground validation data, which is rarely available for old aerial images. As such, our method of aerial image processing could only be used in other similarly documented cases with available old imagery, while more generally applicable proxies are needed to expand investigations over larger areas.
The main advantage of Landsat satellite data compared to other satellite programs is that it provides the longest time series of comparable data (TM5, ETM + 7, OLI8 sensors) from 1984 to the present (Zhu 2017). Landsat resolution (30 m) does not allow the detection of fine-scale delineations like flark boundaries, but the abundance of finescale features with contrasting spectral properties within Landsat pixels, like flarks and Sphagnum lawns in the case of aapa mires, affects reflectance. In the present study, blue, green, red, NIR, SWIR1, and SWIR2 bands were highly correlated with the flark area. The wet flark surfaces had low reflectance in all spectral regions. However, NIR had the best contrast between the wet flarks and the full cover of Sphagnum moss lawns. In general, individual Landsat bands performed equally well or better than most spectral indices in detecting variations in the flark area from aerial images. A corresponding conclusion was reached by Shibayama and others (2009Shibayama and others ( , 2011Shibayama and others ( , 2012 in monitoring the development of rice plants in paddy fields and Mondejar and Tongco (2019) in examining the identification of water bodies. These studies showed that the NIR wavelength range was supe-Ongoing Fen-Bog Transition rior to spectral indices. However, in our case, the GDVI worked almost as well as NIR.
In our study area, Sphagnum-dominated areas stood out in the landscape with high NIR reflectance. Similar results have been previously reported from mires of Hudson Bay lowlands (Whiting 1994) and northern Finland (Middleton and others 2012), as well as from northern Scottish blanket bogs (Lees and others 2020), while some laboratory-based studies have indicated a lower NIR reflectance of Sphagnum mosses than vascular plants (Bubier and others 1997;Harris and others 2005;Pang and others 2020). In nature, Sphagnum mosses grow in mixed communities with vascular plants, and their combined spectra are also affected by leaf area, leaf angle, and litter (Asner 1998). The vegetation types that expanded in the flarks in our study area were Sphagnum carpet and lawn vegetation. They had high NIR reflectance in comparison with wet flarks and coniferous forests, and to water bodies and wooded mire margins, that is, the basic features of boreal mire landscapes. However, vegetation types and their spectral properties, as related to the investigated case, need to be well recognized, and our findings may not be applicable as such in other peatland types.
The moisture content of Sphagnum mosses affects their NIR reflectance, and water loss generally increases it (Bryant and Baird 2003;Harris and others 2005;Pang and others 2020). Sphagnum mosses have a high water-storage capacity, but as poikilohydric plants, they are susceptible to desiccation during droughts (van Breemen 1995;Harris and Byant 2009). This is especially true for carpet and lawn species of Sphagnum, while hummock species resist drying better (Rydin 1985; Há jek and others 2014). In our case, Sphagnum carpet and lawn vegetation had relatively high NIR reflectance despite a high moisture content. However, it is likely that varying moisture conditions contributed to the fluctuation of the NIR level present in our time series. Other factors that may cause fluctuation include slightly different phases of the growing season and artifacts resulting from atmospheric effects, sensor view angles, and radiometric differences between sensors (Jones and Vaughan 2010; Sulla-Menashe and others 2016). However, compared to visible light, there is less absorption and scattering in the atmosphere in the NIR range, allowing higher dynamic range and better contrasts in the images (Smith 1997;Pinto and others 2020). In our case, the differences in timing of selected Landsat images during the growing season could be shown to have contributed to variations in NIR reflectance (Figure 8). When applying the Landsat NIR band in flark change detection, variation in the reflectance levels must be controlled for comparability of the images and applied with appropriate preprocessing.
Despite the fluctuation, our Landsat NIR time series showed a clear and consistent pattern of change that was rarely influenced by potential artifacts. Through the time series, the Sphagnum bog areas of our study site and the reference open bog had a consistently high NIR reflectance. In the beginning of the time series, the flark fen had low reflectance, similar to a coniferous forest. This setting changed drastically, as the NIR reflectance of flark fen approached the values of the Sphagnum mire and reference raised bog by 2010. In the same timeframe, the NIR reflectance of coniferous forest remained at the same low level, indicating that the observed trend in the aapa mire was not caused by the artifact distortion of data. The shift to high reflectance can be explained by the increase in Sphagnum cover and the corresponding decrease in open water in the flarks, which was verified by vegetation resurvey.
The Landsat NIR time series result suggests that the vegetation change included a nonlinear shift that only took place very recently in the 2000s. However, the aerial image series did indicate a decrease in the flark area already between 1947 and 1988, and more rapidly between 1988 and 1997. It is probable that the Landsat NIR method and the aerial image method indicate different thresholds and processes of flark changes that may have different timings. For example, gradual infilling of flarks by vegetation was likely reflected in our Landsat analysis well before the flarks were classified as infilled in our aerial image analysis, which depended on the threshold of flark identification. Importantly, NIR levels continued to rise after the infilling of flarks, probably indicating the continuing change of vegetation growing higher above the water level. In every case, recognition of flarks is central to definition of aapa mires and the phenomenon of their infilling by Sphagnum mosses defines a major ecosystem-scale change. Methods with different resolutions and approaches of both classification and continuous spectral data can complement the assessment of changes in aapa mires.
In our results, NDVI (and several other vegetation indices) showed disparate correlations between years. In 1988, the correlation between NDVI and flark area was negative. This could be explained by the low NDVI of sparsely vegetated flarks. In 2019, however, the correlation was positive, which could be related to recently increased vegetation cover in the flarks. The non-significant correlation in 1997 would then be explained by the mire being in a transition phase. These altering correlations are affected by differences in the growing season stage and can hardly be interpreted conclusively in our study, as only two time points had ground validation. In every case, our results call for caution in the interpretation of NDVI in peatland monitoring, particularly when open water is involved. In general, spectral indices have been developed to better identify certain factors, such as vegetation greenness or the presence of water, and they also aim to stabilize disturbances caused by the atmosphere, viewing angle, etc. (Jones and Vaughan 2010;Lillesand and others 2015). However, there are cases when single-channel data worked better, as observed in our study, and by Shibayama and others (2009Shibayama and others ( , 2011Shibayama and others ( , 2012 and Mondejar and Tongco (2019), for example, regarding the NIR range of the spectrum.

CONCLUSIONS
In this study, we applied multi-resolution remote sensing time series for detecting changes in vegetation and hydrotopographical patterns in a single aapa mire complex with long-term ground validation data. The flark fen habitat declined drastically in our study area between 1959 and 2018. The open water flarks were infilled by Sphagnum mosses, and the central fen zone became narrower, supporting the interpretation of ongoing fen-bog transition, while the vegetation still indicated weak minerotrophy. The legacy of former flark pools was reflected in the modern pattern of Sphagnum vegetation types, as evident in the classified high-resolution UAV image. The observed changes conform to the expected trajectory of autogenic peatland succession, while recent warming and catchment disturbance may have accelerated the process. This study adds another case among several recent studies reporting progressive changes in northern fens that may lead to the development of new bogs in future and potentially increase the importance of northern mires as global carbon sinks.
Our results showed that the Landsat NIR band had especially high performance with solid theoretical grounds in detecting changes between the main vegetation types in aapa mires. Thus, the trend in NIR reflectance may reveal the process of the fen to bog transition. NIR is suggested for monitoring aapa mire complexes, and the method may be applicable to other mire types with clearly contrasting patterns between vegetative and water surfaces. Although high-resolution aerial image interpretations are unchallenged in the precise localization of mire patterning, Landsat images can provide relevant indications of changes between mire type zones in mire complexes. However, further development of this method is needed, and we are currently testing it in several other aapa mire complexes, where Sphagnum expansion has occurred according to aerial images (unpublished data). In future, new satellite sensors with higher spatial and spectral resolution (for example, Sentinel or PRISMA) will enable the more accurate mapping of vegetation and hydro-topography, although the Landsat archives will retain their value in long time scale assessments. Simultaneously, UAV applications are providing high accuracy that is approaching the level of ground validation for many applications. Automated and cloud-based methods (Murray and others 2018; DeLancey and others 2019; Mahdianpari and others 2020) are likely to become more common, but they do not reduce the need for detailed case studies, ground validation, and an understanding of ecological realities.

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 Ongoing Fen-Bog Transition 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 h ttp://creativecommons.org/licenses/by/4.0/. DATA AVAILABI LI TY Kolari, T. H. M., Sallinen, A., Wolff, F., Kumpula, T., Tolonen, K., & Tahvanainen, T. (2021). Vegetation and GIS data examined in the paper ''Ongoing fen-bog transition in a boreal aapa mire inferred from repeated field sampling, aerial images, and Landsat data'' [Data set]. Zenodo. https://d oi.org/10.5281/zenodo.4889010.