Biogeophysical impacts of peatland forestation on regional climate changes in Finland

. Land cover changes can impact the climate by in-ﬂuencing the surface energy and water balance. Naturally treeless or sparsely treed peatlands were extensively drained to stimulate forest growth in Finland over the second half of 20th century. The aim of this study is to investigate the bio-geophysical effects of peatland forestation on regional climate in Finland. Two sets of 18-year climate simulations were done with the regional climate model REMO by us-ing land cover data based on pre-drainage (1920s) and post-drainage (2000s) Finnish national forest inventories. In the most intensive peatland forestation area, located in the middle west of Finland, the results show a warming in April of up to 0.43 K in monthly-averaged daily mean 2 m air temperature, whereas a slight cooling from May to October of less than 0.1 K in general is found. Consequently, snow clearance days over that area are advanced up to 5 days in the mean of 15 years. No clear signal is found for precipitation. Through analysing the simulated temperature and energy balance terms, as well as snow depth over ﬁve selected subregions, a positive feedback induced by peatland foresta-tion is found between decreased surface albedo and increased surface air temperature in the snow-melting period. Our modelled results show good qualitative agreements with the observational data. In general, decreased surface albedo in the snow-melting period and increased evapotranspiration in the growing period are the most important biogeophysical aspects induced by peatland forestation that cause changes in climate. The results from this study can be further integrally analysed with biogeochemical effects of peatland forestation to provide background information for adapting future forest


Introduction
Climate response to anthropogenic land cover change happens more locally and occurs on a much shorter time scale compared to global warming due to increased greenhouse gases (GHG) (IPCC, 2013).The influences on the climate from the biogeophysical effects caused by land cover changes can enhance or reduce the projected climate change (Bathiany et al., 2010;Bonan, 2008;Feddema et al., 2005;Gálos et al., 2011;Göttel et al., 2008;Ge and Zou, 2013;Pielke et al., 2011Pielke et al., , 1998;;Pitman, 2003).Especially for the climate impacts of past large-scale afforestation, studies show that the most obvious effects of the increase of forests in boreal areas are warming during snow-cover periods due to decreased surface albedo and cooling in summertime from increased evapotranspiration (ET) in tropical areas with sufficient soil moisture (Bala et al., 2007;Betts, 2000;Betts et al., 2007).
Vast areas of naturally treeless or sparsely treed peatlands have been drained to grow forests for timber production in northern European countries (Päivänen and Hånell, 2012).In Finland, it is the dominant land cover change over the last half century due to the high fraction of pristine peatland and the need for timber production.The total peatland area of Finland was estimated to Published by Copernicus Publications on behalf of the European Geosciences Union.be 9.7 million ha in the 1950s (Ilvessalo, 1956).In the beginning of the 2000s, the area of drained peatland for forestry was estimated to be 5.7 million ha by Minkkinen et al. (2002) and 5.5 million ha by Tomppo et al. (2011).The area of drained peatlands is unlikely to increase further because no more public subsidisation is given for the first-time drainage of peatlands and the increased awareness of natural conservation (Metsätalouden kehittämiskeskus Tapio, 1997).The area of restored mires was 15 000 ha between 1990 and 2008 (http://www.biodiversity.fi/en/indicators/mires/mi17-mire-restoration) (Kaakinen and Salminen, 2006).However, land cover change is not only a result of human land-use activities but can also be a consequence of climate change.Global warming in the future is also considered to be a factor that affects boreal peatland through water-level drawdown due to increased ET (Laiho et al., 2003;Laine et al., 1995).
Attention has been paid to the climate effects of peatland forestation.A decrease in the local night-time minimum temperature during the growing season was observed roughly for the first 15 years after drainage (Solantie, 1994).The reason for this nocturnal cooling phenomenon is the insulation of lower soil layers from the atmosphere by dry peat.Therefore, the heat flux from drained peat soil can not compensate for the radiative cooling at the surface, which leads to a drop in daily minimum temperature (Venäläinen et al., 1999).On a longer time scale, the growing forest on formerly open peatlands leads to a decrease in surface albedo.The reasons for this are the darker tree cover in comparison to the lighter moss/grass cover in the snow-free period and the partial snow cover in forest areas compared to the full snow cover in open areas in the snow-cover period.This increases the daily maximum temperature due to an increase in the absorption of short-wave radiation (Solantie, 1994).Consistent results on the seasonal cycles of surface albedo and net surface solar radiation due to peatland forestation were found by Lohila et al. (2010), based on measurement data at two pairs of drained and undrained peatland sites located in the south and north of Finland.The results showed a notably decreased surface albedo and corresponding increased net surface solar radiation in springtime.Furthermore, Lohila et al. (2010) indicated the local climate impacts of peatland forestation by investigating long-term  spring surface temperature trends over southern (< 65 • N) and northern (> 65 • N) Finland.The largest positive daytime maximum temperature trend of 0.64 K decade −1 happened in April in southern Finland, where a total of 2.7 million ha of peatlands were drained (Hökkä et al., 2002).The night-time minimum temperature trend through the same period was 0.37 K decade −1 .Lohila et al. (2010) attributed the substantially larger increase in the daytime maximum temperature than in the night-time minimum temperature to the change in surface radiative properties after drainage.
However, these studies about the effects of peatland forestation on climate are based on site-level data or observation-based regional data, which can not attribute the climate impacts to different influencing factors.Specifically, they can not distinguish the local biogeophysical effects from the global climate change due to the increase in GHG concentrations.The climate effects of peatland forestation have not been quantified on a regional scale/country level to investigate the biogeophysical effects in particular.Also, the magnitude and pattern of land-use change effects on climate depend on regional conditions such as soil property, topography, etc. Information from regional studies is essential for the development of future strategies for climate mitigation or forest management.Thus, it is necessary to investigate the effects regionally and systematically.
In recent years, regional climate models have become suitable for simulating regional climate in a fine resolution to resolve small-scale atmospheric circulation (Déqué et al., 2005;Jacob et al., 2001Jacob et al., , 2007;;McGregor, 1997).For this, a regional climate model with a realistic land scheme to interpret more detailed land surface information needs to be applied.
In this study, the long-term climate effects caused by peatland forestation are assessed from two sets of 15-year simulation results with the regional climate model REMO, by using the historical (1920s) and present-day (2000s) land cover conditions.The intention of this study is to understand how peatland forestation in Finland influences regional climate conditions through biogeophysical processes.

REMO climate model
The regional climate model REMO is a three-dimensional hydrostatic atmospheric circulation model developed at the Max Planck Institute for Meteorology in Germany (Jacob et al., 2001(Jacob et al., , 2007;;Jacob and Podzun, 1997).Its dynamical core is based on the "Europa-Modell", the former numerical weather prediction model of the German Weather Service (Majewski, 1991).The land surface scheme (LSS) of REMO mainly follows that of the global atmosphere circulation model ECHAM4 (Roeckner et al., 1996) with several physical package updates (details are shown below).The prognostic variables are pressure, temperature, horizontal wind components, specific humidity, cloud liquid water and ice.REMO is driven by large-scale forcing data according to the relaxation scheme (Davies, 1976).The eight outermost grid boxes at each lateral boundary are the sponge zone.
Because land cover is central to this study, a brief introduction of the LSS in REMO is given below.In REMO LSS, the total area of each model grid box is composed of fractions of land (vegetation cover and bare soil), water (ocean surface and inland lake) and sea ice (Semmler et al., 2004).The biogeophysical characteristics of major land cover classes (Olson, 1994a, b) are described by the following surface parameters: background surface albedo (albedo over snowfree land areas), roughness length, fractional green vegetation cover, leaf area index (LAI; one-sided green leaf area per unit ground area), forest ratio (fr; fractional coverage of trees regardless of their photosynthetic activity), soil waterholding capacity (maximum amount of water that plants may extract from the soil before wilting begins) and volumetric wilting point (percentage of moisture in a soil column below which plants start to wilt) (Hagemann, 2002;Hagemann et al., 1999).The land surface parameters are averaged linearly according to fractional coverage of land cover types within a model grid box, except for the roughness length that is averaged logarithmically (Claussen et al., 1994;Hagemann et al., 1999).As LAI, fractional green vegetation cover and background surface albedo strongly depend on the vegetation phenology, they are prescribed with intra-annual cycles by using a monthly varying growth factor that determines the seasonal growth characteristics of the vegetation (Hagemann, 2002;Rechid and Jacob, 2006).The growth factor for latitudes higher than 40 • north or south is derived from a 2 m temperature climatology (Legates and Willmott, 1990); in other latitudes, the fraction of photosynthetically active radiation is used.
The simple bucket scheme (Manabe, 1969) is used for soil hydrology where the partitioning of surface runoff and infiltration follows the Arno scheme (Dümenil and Todini, 1992).The soil temperature profile from the ground surface to around 10 m deep is described by five soil layers with increasing thickness.The heat conductivity and heat capacity, required in the heat conduction equation for calculating the soil temperature, depend on the soil types (Kotlarski, 2007).The distribution of soil types is from the FAO/UNESCO soil map of the world (FAO/ UNESCO, 1971UNESCO, -1981;;Kotlarski, 2007).
The Arno scheme used for the soil hydrology was further improved by considering the high resolution subgridscale heterogeneity of the field capacities within a climate model grid box (Hagemann and Gates, 2003).The resolution of subgrid-scale heterogeneity is set to be 10 times higher than the model resolution when using the default REMO land cover map-Global Land Cover Characteristics Database (GLCCD) (Loveland et al., 2000;US Geological Survey, 2001).The three parameters in the improved Arno scheme account for the shape of the subgrid distribution of soil water capacities (Beta), subgrid minimum (W min ) and maximum (W max ) soil water capacities.Also, the original annual background albedo cycle was modified by using MODIS satellite data between 2001 and 2004 in order to derive more realistic global distributions of pure soil albedo and pure vegetation albedo, which are then used to compute the annual background albedo cycle with monthly varying LAI (Rechid, 2008;Rechid et al., 2009).

The model domain and land cover data sets
Our model domain covers Fennoscandia, a part of Russia and the northern part of central Europe, and it is centred on Finland (Fig. 1).Typical features influencing the climate of this domain include the North Atlantic Ocean and the Baltic Sea that surround the Fennoscandian countries, many inland lakes located in Sweden and Finland and the relatively high Scandinavian mountain range; the rest of the area has a topography lower than 300 m above sea level.
The default land cover map in REMO is the GLCCD.However, its description of the land cover in Finland is unrealistic.For instance, there is no peatland in Finland in the GLCCD, whereas 7.4 % (22 377 km 2 ) of the land is covered by naturally treeless or sparsely treed peatlands according to the 10th Finnish national forest inventory (FNFI10) (Korhonen et al., 2013).The GLCCD was therefore substituted by the more realistic and up-to-date CORINE land cover map (CLC;2006) for the same model domain in Gao et al. (2014), except for the Russian part where the CLC ( 2006) is not available.Unfortunately, land cover maps describing the land cover conditions of Finland before the most intensive period of peatland drainage in the 1960s are quite limited.Nevertheless, the data collected in the 1st Finnish national forest inventory (FNFI1) provide the possibility for tracing back the land cover condition of Finland in the 1920s (Ilvessalo, 1927;Tomppo et al., 2010).Also, the FNFI10, rather than the CLC (2006), is adopted to describe the land cover condition of Finland in the 2000s, with the aim to avoid the uncertainties in comparing land cover maps with different land cover classification methods and different spatial resolutions.The FNFI1 and FNFI10 land cover maps are post-products that were specially prepared for this study from the respective FNFI field measurement data.The detailed description of the procedures for deriving the FNFI1 and FNFI10 land cover maps is shown in Appendix A. The two FNFI land cover maps are in 3 km resolution and include 10 land cover classes following CLC nomenclature.
In the FNFI maps, the land cover class "peat bogs" is defined as naturally treeless peatland and pine mires where the stocking level is low or the mean height of trees is below 5 m at maturity.Therefore, the shifting from peat bogs to forests represents a major land cover change due to peatland forestation.
In addition to regional inspections, five subregions were selected to represent different land cover change conditions between FNFI1 and FNFI10 (Fig. 1), and the changes of fractional coverage of the 10 land cover classes in those five subregions are given in Table 1.This was done to specifically assess the local climate effects of different intensities of peatland forestation.From subregion1 to subregion4 there is a decrease in the reduction of peat bogs.Subregion1 and subregion2 are two peatland forestation areas located in the middle and south of Finland respectively.In subregion1 and subregion2 there were decreases in the fractional coverage of peat bogs of more than 20 %, and the decreases were mainly compensated by coniferous forest.The decrease in the fractional coverage of peat bogs was 2 % less in subregion2 than that in subregion1, but the increase in the fractional coverage of coniferous forest was 5 % higher in subregion2 than that in subregion1.The total increase in the fractional coverage of forest types was about 16 % in both subregion1 and subre-gion2.Subregion3 is located in the east of subregion1.There was a 12 % decrease in the fractional coverage of peat bogs, but instead of an increase of forests, the fractional coverage of transitional woodland/shrub increased by 14.3 %.Sub- region4 is an area where the most intensive anthropogenic activities have occurred in the five subregions.There was a 14 % decrease in the fractional coverage of forest types and a 3.8 % decrease in that of peat bogs, with a 5.7 % increase in the fractional coverage of artificial areas and a 10.5 % increase in that of agriculture areas.Subregion5 is an area with an 8.64 % increase in the fractional coverage of peat bogs and a 16.3 % decrease in the fractional coverage of forest types.Herein one should notice that some uncertainties may arise from sampling in the FNFI1 and FNFI10 data.This applies especially to FNFI1, where the distance between inventory lines was as high as 26 km.Therefore, subregions that are smaller than 100 km × 100 km may not be sufficient to represent the actual land cover changes spatially.However, the dynamics of the local effects of land cover changes on climate can not be detected when averaging climate signals over large areas with diverse land cover changes.Therefore small subregions, which cover a range of land cover change intensities, are chosen to reflect local climate impacts due to different land cover changes.Moreover, the FNFI data only cover the land surface in Finland without considering inland lakes.Therefore, the land-sea mask in the model domain is adopted from the CLC (2006).In addition, the land cover conditions of the area outside Finland in the model domain are the same as those in Gao et al. (2014), i.e. based on the CLC ( 2006) and the GLCCD, and thus identical in both simulations.
In order to make the land surface parameters more suitable for this study, several modifications in REMO LSS were done.Details of those modifications are documented in Appendix B.

Experiment design
Two simulations were conducted with the FNFI1 and FNFI10 land cover maps, representing the land cover conditions before and after peatland forestation activities in Finland respectively.The simulations were driven with 6-hourly lateral boundary conditions from ECWMF ERA-Interim re-analysis data (Simmons et al., 2007) from 1 January 1979 to 31 December 1996.The 18-year forward runs were preceded by 10-year (1 August 1979-1 January 1990) simulations in order to stabilise the deep soil temperatures and soil moistures.The last 15 years (1 December 1981-30 November 1996) out of the 18-year forward simulations were adopted for further analysis.The analysed period starts from 1 December in order to keep all 3 winter months continuous.The simulated first 1.5 years were excluded in order to minimise the influences of the initial boundary conditions on simulated climate conditions, which have a much quicker adaptation speed than deep soil temperature.The model grid is in an 18 km resolution horizontally and extends over 27 vertical levels (up to 25 km).The model time step was set to 90 s and the time steps of output variables are 6-hourly for 3-D variables and hourly for 2-D variables.Daily data covering 24 h are processed from 18:00 UTC on the previous day to 17:00 UTC on the current day.For 6-hourly data, 18:00 UTC on the previous day and 00:00 UTC, 06:00 and 12:00 UTC on the current day were used for daily values.For this study domain, the growing season and the dormancy season cover the period from May to October and from November to April respectively.

Results
The land cover change effects on regional climate conditions in Finland are analysed based on the differences in climate variables between the post-drainage and pre-drainage simulations (FNFI10-FNFI1).This "delta change approach" is adopted to eliminate the uncertainties related to model bias (Gálos et al., 2011;Jacob et al., 2008).

Effects on climate over Finland
The differences in monthly-averaged daily mean 2 m air temperature (T 2 m ) are quite heterogeneous temporally and spatially.T 2 m differences are most prominent in springtime and summertime (Fig. 3).The most noticeable difference in T 2 m , up to 0.43 K, takes place in the most intensive peatland forestation area in the middle west of Finland in April.The warming is also evident in February and March, with differences of 0.2 K in this area.However, T 2 m turns to show a slight cooling, generally less than 0.1 K, in a few parts of this area from May to October.There are also two regions in northern Finland that show opposite changes compared to the peatland forestation area in the middle west of Finland with cooling in the spring and warming in the growing season.This is because of decreased forest cover and increased fraction of peat bogs in those two areas from FNFI1to FNFI10-based land cover maps.An increase of less than 0.2 K is seen in T 2 m in the southeast of Finland in July and August as well as in the very south of Finland throughout the growing season, which is mainly due to the change from mixed forest to coniferous forest and the increased artificial areas respectively.The 15-year averaged monthly precipitation shows only small differences, less than 10 mm month −1 , in varied patterns in the model domain from April to August (not shown).
The snow clearance day is also an important indicator of springtime climate change at high latitudes (Peng et al., 2013).Therefore, the snow clearance day for each grid box in Finland is determined for the 15 years.The snow clearance day is defined here as the first day after which the total number of snow-covered days does not exceed the total number of snow-free days, and the selection of this day ends before midsummer in a year.The differences between the 15-year averaged snow clearance days of the two simulations (Fig. 4) show almost the same pattern as the differences in T 2 m in April (Fig. 3).In the peatland forestation area in the middle west of Finland, the snow clearance days are mostly advanced by 0.5 to 3 days and, in a few grid boxes, advanced by up to 5 days in the 15-year mean.The two small areas in the north of Finland with reverse land cover changes in comparison to peatland forestation show up to 2-day delays in general.In the very south of Finland, the snow clearance days are also generally advanced in accordance with the warming seen in T 2 m , but delayed in several scattered grid boxes due to increased fraction of artificial areas at the expense of forests.

Effects on climate over five subregions
T 2 m , precipitation and several closely related climate variables (surface albedo, net surface solar radiation, snow depth, ET) for the five subregions were processed into 11-day running means to reduce the influence of day-to-day variations.The differences between the simulations in each of the regionally averaged climate variables were further averaged over the 15 years (Fig. 5).The date information herein (day of year, DOY) represents the middle contributing day of the 11-day averaging period.
T 2 m of subregion1 shows a warming of 0.1 to 0.2 K from February until the end of March and an evident peak of increase from early April to early May (from DOY 95 to DOY 125) that reaches a maximum of 0.5 K in late April.T 2 m of subregion2 has the same development as subregion1 throughout the whole year, but the warming is much smaller and the biggest difference, only 0.12 K, occurs in the beginning of April.This is consistent with the differences in snow depth.The snow-cover period in subregion2 is shorter along with an earlier maximum difference in snow depth.Moreover, those characteristics of the differences in snow depths are in qualitative agreement with the differences in surface albedo because snow is the key factor that controls the surface albedo in the snow-cover period.From the beginning of May to the beginning of October, T 2 m shows a cooling of less than 0.1 K in subregion1 and subregion2 because the cooling caused by ET exceeds the warming caused by the slightly lower albedo.The variability of the differences in net surface solar radiation in the growing season is induced by the variability of cloud cover rather than surface albedo.In November, December and January, the differences in T 2 m vary in both directions.At high latitudes, incoming solar radiation is quite small and cloud cover fraction is high in late autumn and winter.Therefore, the differences in surface albedo are not able to induce differences in net surface solar radiation in this period.Instead, the surface air temperature is sensitive to changes in the long-wave radiation balance that may lead to atmospheric air temperature inversion under a clear sky, manifesting itself as extreme cold surface air temperature.Thus, the variability of the differences in cloud cover caused by short-term variations in the climate contributes to the varied differences in T 2 m in this period.

Discussion
The differences in T 2 m for subregion3 show a warming of less than 0.1 K from DOY 91 to DOY 120 but also a warming in an even smaller magnitude throughout the growing season.The difference in surface albedo in subregion3 is close to 0, although the difference in snow depth is similar to that of subregion2 but with a time lag of around 15 days in the most intensive point.In subregion4, the snow depth shows a quite small increase from the beginning of January until the end of March, which is consistent with the increase in surface albedo and explains the slight decrease of up to 0.1 K in T 2 m from the middle of February until the end of March.Subre-gion5 displays the opposite characteristics compared to sub-region1 and subregion2 for all the investigated variables.The absolute differences in snow depth of subregion5 are smaller than those of subregion1 but larger than those of subregion2.Because subregion5 is located in the north of Finland, the biggest difference in snow depth occurs later than that of sub-region1.The magnitude of the maximum differences in T 2 m in the snow-cover period of subregion5 also lies between that of subregion1 and subregion2 and happens later than that of subregion1.
The differences in T 2 m in the growing season depend on the surplus of energy balance terms where ET manifests itself as latent heat flux.In general, the increase of ET in sub-region2 is slightly higher than that in subregion1.As a consequence, the decrease of T 2 m in subregion2 is slightly larger than that in subregion1 during the growing season when the albedo difference is quite small.The decreased ET and the slightly decreased surface albedo together result in a slight warming during the growing season in the other subregions.The extents of warming in the other subregions follow the magnitudes of the decreased ET because the differences in surface albedo are almost the same in the growing season.
Precipitation has higher variability than ET throughout the year in the five subregions.In general, the differences in precipitation are much larger in the growing season than in the dormancy season, when they are close to 0 mm day −1 .In the growing season, the increase in precipitation of subregion1 occurs during a longer period and has a larger magnitude than that of subregion2.There are slight increases in the precipitation in subregion3 and subregion4, whereas the precipitation of subregion5 shows a decreasing tendency in the growing season, with the biggest differences less than 0.2 mm day −1 .
Furthermore, the maximum and minimum differences of grid pointwise and regionally averaged 11-day running mean of T 2 m over 15 years for subregion1 were investigated as complements to the regionally averaged 15-year mean differences (Fig. 6).T 2 m shows a maximum difference in grid pointwise of nearly 2 K in the snow-melting period over the 15 years, which is 1 K higher than the maximum difference in regionally averaged T 2 m over the 15 years and 4 times that of the 15-year mean of regionally averaged T 2 m .The timings of the three kinds of maximum differences in spring deviate from each other by 3 to 10 days.The minimum differences show only a small deviation between the grid pointwise and regional mean values over the 15 years.During the snowmelting period, the minimum differences of regionally averaged T 2 m is above 0, but not that of the grid pointwise T 2 m .The springtime differences between regional mean and grid pointwise extremes elucidate that, even within one subregion with homogenous characteristics related to peatland forestation, the spring warming of T 2 m is temporally and spatially heterogeneous.This implies that local effects are more pronounced than the regional and temporal statistics can reveal.
For the rest of the year, the differences between the maximum (minimum) of the grid pointwise and regionally averaged T 2 m are small and of a more regional nature.In the period between November and January, the large variations of maximum (minimum) T 2 m are contributed by the inversion effects due to short-term variations in the climate.
Additionally, for a more thorough understanding of the relationships between spring warming and albedo changes in the snow-cover period due to peatland forestation, two correlation relationships were investigated over the 15 years for subregion1 (Fig. 7).One is between the maximum temperature difference day (DOY) and the maximum surface albedo difference day (DOY).The other is between the inflection day of total albedo (the day when surface albedo just finishes a fast decrease from its wintertime level; DOY) and the snow clearance day (DOY).The maximum temperature difference days match to maximum albedo difference days in 6 years, and the rest of the years generally show a delayed maximum temperature difference day compared to the maximum albedo difference day, with a maximum deviation of 14 days.In general, the snow clearance day correlates well with the inflection point of surface albedo.For most years, the differences are less than 6 days, but 3 years show differences up to around 20 days.In those years, sporadic snowfall with a small accumulated snow depth cannot really introduce differences in total surface albedo over the subregion but influences the determination of the snow clearance day.

Relationships between the changes in biogeophysical aspects and the impacts on climate
To assess the generality of the causal relationships between land cover changes and climate variables, the spatial correlations between changes in the two surface energy balance relevant variables, surface albedo and ET, and T 2 m are investigated.Consequently, the spatial correlations between changes in surface albedo and ET and changes in the surface parameter values are also explored.The correlations with fractional green vegetation cover is not shown in Fig. 8 because LAI and green vegetation ratio are both modulated with the monthly varying growth factor by the same scheme, and they are highly correlated (Pearson correlation coefficient, r 2 = 0.984 for March, r 2 = 0.674 for June).Monthly means of 15-year averaged changes in March and June are selected to represent springtime and summertime respectively.The changes in T 2 m are in accordance with the changes in surface albedo in March (Fig. 8a), which is almost linearly correlated with the changes in LAI (Fig. 8c) and forest ratio (Fig. 8e).The changes in T 2 m are linearly correlated with the changes in ET over most of the area in June (Fig. 8b).In general, the changes in ET are also correlated with the changes in LAI (Fig. 8d), roughness length (Fig. 8f) and forest ratio (yearly constant, not shown), despite the influences from droughts that may happen in late summer.Overall, the changes in surface albedo and ET are closely dependent on the changes in land surface parameters, which are induced by the changes in fractional coverages of land cover types in the five subregions (Table 1).The changes in T 2 m are mainly modulated by the changes in surface albedo and ET in spring and summer respectively.Some grid boxes located in the southeast of Finland, where mixed forest was substituted by mainly coniferous forest, show deviations in the correlations with LAI (marked by yellow circles in Fig. 8b, c, d).In this area, LAI increased with almost no change in forest ratio, which led to a relatively smaller decrease in surface albedo compared to other areas with the same magnitude of changes in LAI in March; the ET-induced cooling is outweighed by the albedo-induced warming, which causes a slight warming in June.In the following summer months, July and August, the ET-induced cooling typically gets smaller because of surface water limitation and consequent warming.

Biogeophysical impacts of peatland forestation on regional climate
Surface albedo shows a notable decrease in peatland forestation areas during the snow-cover period and a slight decrease in the growing season, whereas LAI, roughness length, fractional green vegetation cover and forest ratio increase throughout the year after peatland forestation.Those changes lead to an increase in springtime T 2 m , which occurs locally in accordance with the decrease in surface albedo.In the growing season, an increase in ET related to the increased LAI and fractional green vegetation cover leads to more energy consumed by latent heat flux than gained by slightly lower albedo.Additionally, higher roughness length can play a role by increasing turbulent mixing and consequently the magnitudes of turbulent fluxes.Thus, the scattered differences in precipitation in summer are contributed to more convective structures, while for the rest of the year the precipitation is basically controlled by large-scale meteorology.From the analysis of the results in the five subregions, the differences in the climate variables show that their magnitudes depend on the extent of land cover changes, while the timings of the extremes mostly depend on geographical locations (latitudes) that define the radiation balance through the seasonal cycle.
Results also illustrate a positive feedback induced by peatland forestation between lower surface albedo and warmer T 2 m in the snow-melting period.The warming caused by lower surface albedo in the snow-cover period due to more forest leads to a quicker and earlier snow melting; meanwhile, the surface albedo is reduced and consequently the surface air temperature is increased.Additionally, the maxi- mum difference in the grid pointwise 11-day running mean of T 2 m in the spring warming period over the 15 years reaches 2 K in subregion1, which is 4 times 15-year mean of the corresponding regionally averaged values.This illustrates that the spring warming effect from peatland forestation is highly spatially and temporally heterogeneous.

Comparison with observation-based results
To examine the realism of the simulated effects on surface temperature in springtime from peatland forestation, linear temperature trends over 40 years  were calculated for March and April based on monthly mean daily maximum (T 2 m,max ) and monthly mean daily minimum (T 2 m,min ) surface temperatures over Finland from an E-OBS gridded observational data set in 0.25 degree resolution (Haylock et al., 2008) (Fig. 9).The significance of the trends was tested with the Student t test.Both T 2 m,max and T 2 m,min have increased in March and April over the 40 years, and the increases of temperatures in March are stronger than those in April.The major areas of peatland forestation, sub-region1 and subregion2, are highlighted and statistically significant (p < 0.1) in the trends of T 2 m,max in both March and April but not shown in the trends of T 2 m,min .In springtime, the trend of T 2 m,max is influenced by the albedo-induced temperature changes locally, while the trend of T 2 m,min is more related to the general climate change caused by global GHG increases.Thus, the local effects in the trends of T 2 m,max suggest that our modelled results show qualitatively a good correspondence to observational data.However, it is difficult to compare the exact magnitudes and patterns of temperature changes because observational data contain contributions from other factors: for instance, the effects of climatic teleconnections from land cover changes in surrounding areas of Finland and short-lived climate forces such as aerosols and reactive trace gases (Pitman et al., 2009).Furthermore, regional averaged difference in the simulated 11-day running mean net surface solar radiation of subre-gion1 (Fig. 5d) agrees well with the observed differences in daily mean  net surface solar radiation (Fig. 4 in Lohila et al., 2010) between open peatland and forest sites located in southern and northern Finland.The maximum differences in the observed net surface solar radiation at nutrientrich sites are 40-45 W m −2 (on DOY 70) in the south and 80-90 W m −2 (on DOY 110) in the north of Finland.At nutrientpoor sites, the maximum differences are 30-40 W m −2 (on DOY 80) in the south and 60-70 W m −2 (on DOY 115-120) in the north of Finland.The maximum difference in the simulated 11-day running mean net surface solar radiation averaged over subregion1 is 6.5 W m −2 (on DOY 107).The timing of the maximum difference in our simulated results for subregion1 falls within the range of that in the observed data.The much smaller magnitude of the maximum difference in the simulated results could be explained by the fact that only around 20 % of the land was transformed from peat- land to forests in subregion1.The maximum difference in net surface solar radiation is caused by the advanced snow clearance day due to peatland forestation.The differences in surface albedo is biggest between snow-covered peatland surface and non-snow-covered forest surface, i.e. the maximum difference of surface albedo is mostly dependent on snow albedo.Snow albedo has a negative linear correlation with forest ratio (Fig. B1).Assuming that the entire land of sub-region1 would have been changed from peatland to forests, the maximum difference in net surface solar radiation could be estimated to be 5 times larger, i.e. 32.5 W m −2 , which is within the range of observations.Moreover, the evolution of the differences in both simulated and observed net surface solar radiation in spring can be divided into three phases: a slow increase, a quick increase and a quick drop.For the simulated net surface solar radiation, the slow increase occurs from the beginning of January until the end of March and appears to be mostly induced by the differences in snow depth of the land cover classes.The following quick increase occurs in a much shorter period in April, within 10 to 20 days.The quick drop of the differences in net surface solar radiation follows the strong decrease of snow cover.The quick increase and quick drop are mainly attributed to snow melting, which is very sensitive to warmed air temperature.After that, the difference of net surface solar radiation is much smaller in the growing season.Strong variability of the net surface solar radiation differences in summertime can be attributed to the relatively short averaging time (15 years) that does not even out the impact of periods of high and low cloudiness.

Perspectives to improve land-atmosphere interactions modelling
Our results show that local climate changes due to peatland forestation in Finland can be mainly attributed to the impacts of changed surface albedo and ET on surface temperature, whereas no strong influences on precipitation is found.Future studies for improving understanding of biogeophysical impacts on regional climate of peatland forestation and other land cover or land use intensity changes could focus on the following issues: the parameterisation of albedo, the soil hydrology scheme, and the implementation and accuracy of land cover maps.
Although the maximum background albedo values of FNFI land cover classes in this study are broadly consistent with the summertime albedo values derived especially for two observation stations in Finland in Kuusinen et al. (2013), the estimated albedo for land cover classes in high-latitude areas show variations in a range of studies.The mean summertime albedo for coniferous forest is only 0.079 in Hollinger et al. (2010), while it is 0.119 in our study.We used a summertime albedo for broad-leaved forest of 0.146, which is higher than the albedo values for deciduous in Kuusinen et al. (2013) but still lower than 0.156 for aspen in Betts and Ball (1997) and 0.152 for deciduous in Hollinger et al. (2010).The cropland albedo is 0.189 in Hollinger et al. (2010), much higher than the cropland albedo of 0.156 used in our study.In the middle boreal zone of Finland, the albedo of peat bogs and the albedo of forest are, on average, 0.145 and 0.115 in Solantie (1988) respectively.Thus, compared to those values, our lower albedo for peat bogs and higher albedo for forest (even only considering coniferous forest) may underestimate the warming effect contributed by more absorbed solar radiation in the non-snow-covered period.However, it is hard to estimate the overall influence on surface temperature because ET may be enhanced from increased net surface solar radiation.Furthermore, even albedo values of the same land cover class could be different in different parts of Finland. In Solantie (1988), the mean albedo of barren bogs in southern Finland and of the concentric raised bogs in the middle of Finland is only 0.128.Recent studies show that forest albedo is influenced by stand density and understory in different sites (Bernier et al., 2011).
In wintertime, the snow albedo scheme is much more important than the background albedo in determining the sur-face albedo for high latitudes.The snow albedo scheme in REMO does not adequately represent the complex conditions over forests with the linear dependence on snow surface temperature.Snow properties and canopy conditions, such as snow water content, grain size and snow pack thickness, as well as impurities on the snow surface, have a strong influence on snow albedo (Wiscombe and Warren, 1980).Moreover, there is no vertical structure of forests in REMO where the process of snow intercepted by canopy is crucial (Roesch et al., 2001).The canopy of forests is also important in causing a night-time warming by the shelter effect in areas with successful peatland forestation after about 15 years (Venäläinen et al., 1999).
Besides, soil moisture affects ET and precipitation, thus playing a vital role in energy partition and influencing surface temperature (Hagemann et al., 2013).The simple bucket soil moisture scheme used in this study is insufficient to represent the complex soil hydrological processes (Hagemann and Stacke, 2014).Also, the subgrid variability of soil saturation within a model grid box is taken into account as one-third the model resolution in the simple bucket hydrology scheme in REMO LSS for this study, which is restricted by the 3 km resolution of the FNFI land cover maps.This can lead to underestimation of the surface runoff because the differences between the two surface parameters, W max and W min , are smaller over the model domain compared to those when using a 10 times finer resolution to represent the subgrid hydrologic heterogeneity with the GLCCD or CLC (2006).
Furthermore, land surface parameters are allocated according to distributions of land cover types in land surface scheme.Spatially more explicit land cover maps with a parameter set tailored for the study area could reduce the uncertainties in simulation results of climate models from the source.

Biogeochemical aspects related to peatland forestation
Peatland is a significant source of CH 4 , and the CH 4 emission rate is sensitive to temperature, water table level, plant root depth, soil nutrition level, etc. (Melton et al., 2013;Turetsky et al., 2014;Lohila et al., 2010).After peatland forestation, the soil water table level goes down, leading to increased CO 2 release at the expense of CH 4 release (Minkkinen and Laine, 2006).As time goes by, carbon sequestration by the tree growth and the formation of a new litter layer could compensate for the carbon loss from peatland.Lohila et al. (2010) combined the radiative forcing effects from the differences of albedo and GHG fluxes due to peatland forestation at sitelevel and showed net cooling at two soil-nutrient-rich sites in the south and north and one soil-nutrient-poor site in the south of Finland.Accounting for such local impacts in a regional climate model requires very sophisticated process descriptions and detailed parameterisation of soil properties.

Y. Gao et al.: Biogeophysical impacts of peatland forestation on regional climate changes in Finland
The biogeophysical impacts of vegetation-climate feedbacks on climate are modest in comparison to the effects of increased GHGs for Europe, but local, regional and seasonal effects can be significant (Wramneby et al., 2010).However, studies with dynamic vegetation models under climate projections with increased GHGs indicate that more carbon will be gained in terrestrial ecosystems at high latitudes by the end of this century (Fallon et al., 2012;Zhang et al., 2014).This is due to an increase in woody plants that induce biogeophysical feedbacks with an earlier onset of the growing season.

Summary
To get a clear picture of the peatland forestation effects on the climate in Finland it is important for future forest management to consider economic aspects and global warming mitigation.In this paper we investigated the long-term biogeophysical effects of peatland forestation on near-surface climate conditions in Finland by using a historical (1920s) and a present-day (2000s) land-use map based on Finnish national forest inventory data in the regional climate model REMO.The differences between the two simulations in surface air temperature and precipitation were examined.The results show that peatland forestation induces a spring warming effect and a slight cooling effect in the growing season, but a varied pattern with less than 10 mm month −1 differences in precipitation over Finland from April to September.The temperature response in spring in simulation results is well in line with that seen in observational maps.In the most intensive peatland forestation area in the middle west of Finland, the monthly-averaged daily mean surface air temperature shows a warming effect of around 0.2 K in February and March and up to 0.43 K in April, whereas a cooling effect of, in general, less than 0.1 K is found from May until October.Consequently, the snow clearance days in model grid boxes over that area are advanced up to 5 days in the mean of 15 years.Furthermore, a more detailed analysis was conducted on five subregions with decreased fractions of transformation from peatland to other land cover classes.The 11day running means of simulated temperature, surface albedo, net surface solar radiation and snow depth, as well as precipitation and ET, were averaged over 15 years.Results show a positive feedback induced by peatland forestation between decreased surface albedo and increased surface air temperature in the snow-melting period.Overall, decreased albedo in the snow-melting period and increased ET in the growing period as a result of peatland forestation are the most important biogeophysical aspects that cause changes in surface air temperature.The extent of these climate effects depends on the intensity and geological locations of peatland forestation.
In the future, with the aim of getting a more precise assessment of the biogeophysical impacts of peatland forestation on regional climate conditions, more accurate land cover maps and land surface parameters are essential.Also, a more robust land surface scheme could enhance the representation of interactions between land surface and climate.

Appendix A: Methods in deriving FNFI land cover maps
The sample of FNFI1 (1921)(1922)(1923)(1924) consisted of inventory lines oriented from southwest to northeast at a distance of 26 km across most parts of the country.The total length of measured lines was 13 348 km, and the total number of assessed land figures was 93 922.In the CLC-classification method, mean tree height and crown cover are two important criteria for classifying land-use classes.However, because crown cover was not measured in FNFI1, the growing stock volume corresponding to crown cover thresholds was estimated using naturally regenerated forests and unditched pine mires in FNFI9 (1996-2003) and in FNFI10 (2004-2010), according to vegetation zone, site type, mean height and dominant tree species.Afterwards, fractions of the 10 land cover classes that were used in this study were derived for the FNFI1 sample in FNFI1 by considering land-use class, estimated growing stock volume classes, mean height, vegetation zone, site type and tree species composition.
For the interpolation, the FNFI1 sample lines were split into slices with 1 km intervals in a S-N direction.The fractions of the 10 land cover classes in each slice on inventory line (1380 m on average) were then used in calculating sample variograms.These sample variograms are fitted into a variogram model to derive kriging predictions using the R version 2.15.2 package gstat (Pebesma, 2004;R Development Core Team, 2011).The block kriging was carried out separately for the fraction of each of the 10 land cover classes with isotropic exponential (or spherical) variogram model and a block size of 50 km × 50 km.A raster map in 3 km resolution was then produced for the coverage of the 10 land cover classes.
In FNFI10, a systematic cluster sample (more details can be found at http://www.metla.fi/ohjelma/vmi/vmi10-otanta-en.htm) of 69 388 plots was measured (Korhonen et al., 2013).The distance between clusters of plots (10-14 plots/cluster) varied between 5 km (in southern Finland) and 11 km (in northern Finland).The classification of the FNFI10 data set was processed in a similar way to the FNFI1 data, with the exception that crown cover thresholds for classifying land-use classes can be used directly in FNFI10 because it is assessed.To derive the 3 km × 3 km grid map, the cluster means of the proportions of the 10 land cover classes were first calculated, and then the same interpolation method was used as for FNFI1.

Appendix B: Modifications in REMO LSS in this study
In order to allocate the surface parameters to appropriate land cover classes, the standard GLCCD land cover classes are related to the 10 land cover classes in the FNFI maps through comparing the definitions of land cover classes (Table B1).Most of the surface parameters follow the built-in parameter values.However, large deviations were found when comparing the parameterised albedo with the observed albedo.Moreover, the method for background albedo parameterisation is not suitable for land-use change studies because the vegetation albedo and the soil albedo maps are both derived from satellite albedo data that were measured in 2001-2004 with respect to land cover over that period.A new method, Land Use Character Shifts (LUCHS), has been proposed for land cover change studies (Preuschmann, 2012).It derives the annual background albedo cycle for certain landuse types in one region from good quality remote sensing data sets -a surface albedo data set and a land cover mask -that are produced in the same time period.Unfortunately, LUCHS is not feasible for high-latitude areas where snow cover prevents the possibility of deriving background albedo values from satellite albedo data.Hence, a simplified method is developed in this study to derive the background albedo values of the 10 land cover classes in FNFI land cover maps.It is based on the assumption that the vegetation albedo map and the soil albedo map in current REMO LSS are feasible to describe the albedo values of the land cover condition in FNFI10 because the two data sets are overlapping in time.• 52 E) Finnish observation station.The station values are estimated by a linear unmixing approach with the land-use and forestry maps in combination with the MODIS BRDF/albedo product (Kuusinen et al., 2013).The derived and observed albedo values show good agreement for peat bogs, mixed forest, transitional woodland/shrub and agricultural areas, as well as for artificial areas.Although the maximum albedo values of coniferous forest and broad-leaved forest in this study are around 0.01 higher than those in Kuusinen et al. (2013), they are reasonable for considering albedo differences between land cover classes.The three land cover classes (natural grasslands, moors and heathland, open spaces) are not found at the two stations; however, they take up only small proportions of the FNFI land cover maps.The snow albedo scheme for calculating the surface albedo during the snow-cover period was also found to require some improvements.When there is snow on the ground, the surface albedo in REMO LSS is a function of background albedo, snow albedo and snow depth.The snow albedo depends linearly on snow surface temperature and fr (Kotlarski, 2007).Based on previous studies (Køltzow, 2007;Räisänen et al., 2014;Roesch et al., 2001), the minimum snow albedo at snow surface temperature T = 0 • C (a min ) and the maximum snow albedo at snow surface temperature T ≤ −10 • C (a max ) of non-forested area (fr = 0) in this study were increased from 0.4 to 0.56 and decreased from 0.8 to 0.68 respectively; in addition, the a min and a max of fully forested area (fr = 1) were both decreased to 0.25 (Fig. B1).For descriptions of the dynamics of snow cover, the interested reader is referred to Kotlarski (2007).
Moreover, the three parameters for describing the subgrid heterogeneity of soil hydrology (Hagemann and Gates, 2003), Beta, W min and W max , were calculated in a subgrid scale of 6 km resolution.It is one-third of the 18 km REMO resolution.The reason for this is that the spatial resolution of the FNFI land cover maps is 3 times lower than that of the default GLCCD land cover map.
Corrections were also made to some of the surface parameters of coniferous forest and mixed forest to obtain a better mutual consistency of the surface parameters for the three forest types.For coniferous forest, the fractional green vegetation cover in the dormancy season, the growing seasons and the forest ratio were set to 0.91, 0.91 and 0.8 respectively, as proposed for Fennoscandia by Claussen et al. (1994).For mixed forest, the fractional green vegetation cover and LAI in the dormancy season were revised to be half of those parameters in the growing season.

Figure 1 .
Figure 1.Orography of the model domain and the five selected subregions (subregion1 -blue; subregion2 -red; subregion3 -purple; subregion4 -green; subregion5 -orange).The inner black frame shows the extent of the relaxation zone from the outer boundary, i.e. the eight outer-most grid boxes in each direction of the model domain.

Figure 2 .
Figure 2. Changes of fractional coverage of the 10 land cover classes in Finland from the 1920s to the 2000s (FNFI10-FNFI1).

Figure 3 .
Figure 3.The 15-year averaged differences (FNFI10-FNFI1) in monthly-averaged daily mean 2 m air temperature in spring and summer months.

Figure 4 .
Figure 4.The 15-year averaged differences (FNFI10-FNFI1) in the snow clearance days over model grid boxes in Finland.

Figure 5 .
Figure 5.The 15-year averaged regional mean differences (FNFI10 -FNFI1) in 11-day running mean of daily mean (a) 2 m air temperature, (b) snow depth (presented as equivalent water), (c) surface albedo, (d) net surface solar radiation, (e) ET and (f) precipitation of the five subregions.

Y.Figure 6 .
Figure 6.Maximum, minimum and mean differences of grid pointwise and regionally averaged 11-day running mean of daily mean 2 m air temperature over 15 years in subregion1.

Figure 7 .
Figure 7. (a) Correlation between maximum temperature change day albedo change day (DOY); (b) correlation between inflection day of to surface albedo just finishes a fast decrease from its wintertime level; DO day (DOY).The plots show regional means over subregion1 for all 15 yea

Figure 8 .
Figure 8. Spatial correlations between (a) changes in monthlyaveraged daily mean 2 m air temperature (T 2 m ) and changes in albedo for March, (b) changes in T 2 m and changes in ET for June and relationships between changes in land surface parameters in REMO LSS following land cover changes and changes in albedo (c, e) (changes in ET, d, f) in the corresponding month.The changes in the grid boxes in selected subregions are shown with coloured dots (subregion1 -blue; subregion2 -red; subregion3 -purple; subre-gion4 -green; subregion5 -orange).The grid boxes in yellow circles show the changes in the southeast area of Finland.
et al.: Biogeophysical impacts of peatland forestation on regional climate changes in Finland

Figure 9 .
Figure 9. Temperature trends over 40 years (1959-1998) for (a, c) monthly mean daily maximum and (b, d) monthly mean daily minimum surface temperatures of March and April.The areas covered with black dots are statistical significant (p < 0.1).

Figure B1 .
Figure B1.Modified snow albedo values in the snow albedo scheme (modified based on Fig. 3.6 in Kotlarski, 2007).

Table 1 .
Changes of fractional coverage (%) of the 10 land cover classes from the 1920s to the 2000s (FNFI10-FNFI1) in the five subregions.

Table B1 .
Translations between the 10 land cover classes in FNFI maps and the GLCCD land cover classes.

Table B2 .
Kuusinen et al. (2013)d vegetation albedo values with standard deviations for the land cover classes in the FNFI maps and the threshold used for each land cover class; the minimum and maximum background albedo values (with standard deviations) in the yearly cycle calculated based on the derived soil and vegetable albedo values are also shown, as well as the summertime albedo values (i.e.maximum background albedo values) of the corresponding land cover classes observed at two Finnish stations (Hyytiälä, Värriö) inKuusinen et al. (2013).modelgridboxes that satisfy a requirement of 80 % coverage of one land cover class in FNFI10, are averaged to represent the soil and vegetation albedo values of that land cover class.The 80 % threshold was decreased to 50 % for natural grasslands, peat bogs and artificial areas, as none of the model grid boxes have an 80 % coverage of those land cover classes in Finland.The derived albedo values and the standard deviations for each land cover class in FNFI maps are shown in TableB2.The maximum background albedos, calculated based on the derived soil and vegetation albedo for FNFI land cover classes, are then compared with the summertime albedo of similar land cover classes for a southern (Hyytiälä; 61 • 51 N, 24 • 17 E) and a northern (Värriö; 67 • 48 N, 27 in