The ERA5-Land Soil-Temperature Bias in Permafrost Regions

ERA5-Land (ERA5L) is a reanalysis product derived by running the land component of ERA5 at increased resolution. This study evaluates ERA5L soil temperature in permafrost regions based on observations and published permafrost products. We find that ERA5L overestimates soil temperature in northern Canada and Alaska, but underestimates it in mid-low latitudes, leading to an average bias of -0.08 °C. The warm bias of ERA5L soil is stronger in winter than in other seasons. As calculated from its soil temperature, ERA5L overestimates active-layer thickness and underestimates near-surface (< 1.89 5 m) permafrost area. This is thought to be due in part to the shallow soil column and coarse vertical discretization of the landsurface model and to warmer simulated soil. The soil-temperature bias in permafrost regions correlates well with the bias in air temperature and with maximum snow height. Review of the ERA5L snow parameterization and a simulation example both point to a low bias in ERA5L snow density as a possible cause for the warm bias in soil temperature. The apparent disagreement of station-based and areal evaluation techniques highlights challenges in our ability to test permafrost simulation models. 10 While global reanalyses are important drivers for permafrost simulation, we conclude ERA5L soil data is not well suited for informing permafrost research and decision making directly. To address this, future soil-temperature products in reanalyses will require permafrost-specific alterations to their land-surface models.

. ERA5L is forced by the atmospheric analysis of ERA5 and hence the assimilated observations indirectly influence the simulations. It is delivered at the same temporal resolution as ERA5 and with a higher spatial resolution of 0.1°. ERA5L is currently available for period 1981-2018, and will eventually extend back to 1950 and be 60 updated to the present time with little delay. Note that at the time of writing, only ERA5L data after 2001 have been released to the public and so this evaluation is conducted using data between 2001-2018.

Snow scheme
A more realistic representation of snow is used in the ERA5 land surface model compared with its predecessor, ERA-Interim. 65 ERA5L uses HTESSEL which treats snow as a single layer above the soil with independent prognostic temperature, mass, density, and albedo (Orsolini et al., 2019). The description of snow processes in HTESSEL is summarized by Dutra et al. (2010) as: (1) liquid water with phase changes coexists with ice in the snow pack and is diagnosed from its temperature, mass, and density (Appendix B1); (2) snow density changes according to overburden, thermal metamorphism and retained liquid water following Lynch-Stieglitz (1994) (Appendix B1); (3) albedo changes exponentially with snow age and is adjusted by 70 vegetation conditions; (4) snow cover fraction depends on both snow water equivalent (SWE) and density (Appendix B2).

Soil scheme
Soil heat transfer in ERA5L is governed by the Fourier law. The thermal effects associated with latent heat are accounted for by following the method of Rouse (1984). However, soil thermal conductivity depends only on moisture content, and the influence of phase change is not represented. The upper boundary condition is given by a heat flux at the ground surface, derived from a 75 weighted average over eight subgrid fractions (or "tiles"). A zero heat flux is assumed at the lower boundary. The ERA5L soil column is discretized into four layers with node depths (layer boundaries) at 0.07 (0-0.07), 0.21 (0.07-0.28), 0.72 (0.28-1.00), and 1.89 (1.00-2.89) m.

Observations and quality control
Soil temperature time series from 639 sites located in permafrost regions were compiled from a variety of sources (Table 1, 80 Figure A1. See station metadata from supplement). Sites consist of both meteorological stations and boreholes. Of these, there are 56 from the China Meteorological Administration (CMA, Wang et al., 2015), 105 from World Data Centers (WDC) in Russia and Ukraine, 219 from Nordicana D, 95 from the Geophysical Institute, University of Alaska Fairbanks (GI-UAF), 10 from the Tibetan Plateau observatory of plateau scale soil moisture and soil temperature (Tibet-Obs) (Su et al., 2011), 60 from multiscale Soil Moisture and Temperature Monitoring Network in the Central Tibetan Plateau (CTP-SMTMN) (Yang et al.,85 2013), 40 from the Global Terrestrial Network for Permafrost (GTN-P, Biskaborn et al., 2015), 28 from National Park Services (NPS) in Alaska (Wang et al., 2018), 16 from the U.S. Geological Survey (USGS, Urban and Clow, 2017;Wang et al., 2018), 8 from HiWATER (Che et al., 2019), and 2 from Boike et al. (2018Boike et al. ( , 2019. The permafrost zone of each site was determined ERA-Interim air temperature; and (4) the 1-km circumpolar permafrost map (CP map) which is derived from a statistical model (Karjalainen et al., 2019a).
Whereas ERA5L, TTOP and CP maps represent permafrost as a boolean variable (i.e. present or absent according to soil temperature), the IPA map and PZI map represent permafrost using either a categorical variable (e.g., continuous, discontinuous, sporadic, or isolated permafrost) or a continuous index (0.01-1) as a proxy to approximately represent the proportion 105 of an area underlain by permafrost (i.e the permafrost extent). Following Melton et al. (2019), we apply a threshold of 50% (corresponding to the continuous and discontinuous permafrost zones) and 0.5 for the IPA map and the PZI map, respectively, to allow for meaningful comparison with the other maps. Values greater than this are considered to represent permafrost areas.
The mean annual ground temperatures (MAGT) from the TTOP and CP maps were also used to evaluate ERA5L.

120
MAAT bias and maximum snow depth (SD max ) were selected as candidate variables to be assessed as possible predictors of ERA5L soil temperature bias (see Eq. 1). SD max was defined as the median of annual maximum monthly snow depth during the period 2001-2018. The surface offset (SO), which quantifies the influence of surface conditions such as snow and vegetation cover (Smith and Riseborough, 2002), is defined here as the difference between MAAT and MAGT of the uppermost soil layer in ERA5L.

125
ERA5L ALT was derived by linearly interpolating the ERA5L soil temperature-depth profiles. The TTOP and CP map were derived using an equilibrium model, and MAGT is given as an average of the entire period (MAGT avg  in any of the four soil layers has an hourly temperature below 0°C for two consecutive years (ERA5L H ); (2) the MAGT of the fourth soil layer is below 0°C for two consecutive years (ERA5L A ).

Detailed permafrost simulation example
Our results show remarkable bias of ERA5L soil temperature in winter that is thought to correlate with snow depth (Figure 2).
For this reason, the suitability of ERA5L soil temperature and the effect of the snow density bias are further investigated using 140 a site-specific simulation example at a densely instrumented location near Lac de Gras (LdG), N.W.T., Canada ( Figure 1A).
This simulation provides an opportunity to evaluate ERA5L soil temperature under different terrain (e.g. vegetation, soil properties) and snow conditions. We used GEOtop 2.0 (Endrizzi et al., 2014), a process-based numerical model, to simulate snow characteristics and soil temperature for ten terrain types between September 2015 and August 2017 as described in detail by Cao et al. (2019a). Snow compaction due to wind effects is considered in 1-D for all terrain types except for the tall shrub site 145 (Pomeroy et al., 1993). The snow-correction factor (SCF) is used to scale modeled snow mass via precipitation. It is used as a lumped variable for representing precipitation bias in the driving reanalysis as well as differences between terrain types that are caused by preferential accumulation and lateral transport by snow drifting. The ERA5 reanalysis and its ten-member ensemble are used as forcing data for the simulation. Note that the CP map only represents permafrost distribution north of 30°N (Karjalainen et al., 2019a), the TTOP map represents the permafrost distribution within the Northern Hemisphere (Obu et al., 2019), and the others represent the permafrost area north of 60°S. Permafrost area from the literature is given with their definition in this study.   (Table 2). Soil temperature is found to have a warm bias in western Canada and Alaska but a cold bias in mid-low latitudes such as the QTP, leading to a near-zero wBIAS of -0.08°C (Figure 3). Among the 932 MAGTs from 331 ERA5L grid cells, 20.7% have an RMSE less than 1°C, 53.5% have RMSE less than 2°C, and 68.9% have RMSE less than 3°C.

155
The following linear model was used to predict ERA5L soil temperature bias in permafrost regions using MAAT bias and snow depth as predictor variables: where wBIAS M AAT is the weighted bias of MAAT. The model was fit using 239 grid cells and has an R 2 of 0.47. Both predictors were found to be statistically significant with p < 0.01. The result suggests that MAAT and snow depth both influence Model (UKESM) (Chadburn et al., 2015).
Averaged MAGTs from the CP and TTOP map were bilinearly interpolated to the observed sites and compared against the 165 observations of the deepest soil layer. We found that the performance of ERA5L is intermediate between the two maps (Table 2).
Whereas Karjalainen et al. (2019b) found that the predictive accuracy of their statistical model was similar between permafrost and non-permafrost regions, our results show ERA5L and TTOP soil temperature agree less with observations in permafrost regions than in non-peramfrost regions (Table 2, Figure 3). In addition to the worse performance of MAAT in these regions, the result suggests that HTESSEL may be less suitable for soil temperature simulation in areas with more prevalent snow and 170 soil freezing. The large warm bias of ERA5L soil temperature during winter ( Figure 2) further supports this notion.

Active-layer thickness and permafrost distribution
While ERA5L is not capable of representing deep ALT, our results show that even for shallow ALT grids, the mean ERA5L ALT (1.67 m) was more than twice the mean observed ALT (0.82 m) (Figure 4). ERA5L ALT is substantially overestimated for most (72/79) of the grids, with wRMSE values up to 0.98 m. Excluding glaciers, the mean near-surface permafrost area of 175 the Northern Hemisphere was estimated as 6.6±0.6×10 6 km 2 based on hourly soil temperature and 9.9±0.5×10 6 km 2 based on MAGT during 2002-2018 (Table 2, Figure 5). ERA5L underestimates permafrost area compared to previous estimations (e.g., Brown et al., 1997;Gruber, 2012;Obu et al., 2019;Karjalainen et al., 2019b). Near-surface permafrost area of ERA5L as defined in this study decreased at a rate of -0.11 (-0.08) ×10 6 km 2 year -1 based on hourly (annual) mean soil temperature.
This corresponds to a loss of 1.7 (1.4) ×10 6 km 2 of near-surface permafrost area since 2002.

Detailed permafrost simulation example
The detailed example simulation indicates that ERA5L soil temperature has warm bias (from 0.95 to 5.48°C) in all terrain types, whereas GEOtop forced by ERA5 and its ten ensemble members show more reasonable results even when SCF = 1 ( Figure 6). More specifically, ERA5L was only found to be suitable in terrain types with significant snow deposition, (e.g. in snowdrifts, tall shrubs, and sedge fen). For all other terrain types, ERA5L showed a significant warm bias during winter and, 185 consequently, in the annual mean. Although the ERA5L results for SWE agreed with GEOtop when driven with the same data (SCF=1), the mean snow depth was approximately 1.53 times that of GEOtop and the snow density was much lower (Table 4).  ERA5L has a number of advantages for permafrost research; it provides a long historical record (back to 1950, eventually), high spatial resolution, and global coverage. While it could be seen to provide an opportunity to study long-term changes of  ERA5L are generally worse in the shallow soil layers (Figure 2). ERA5L does not reproduce ALT well (Figure 4), likely due to its shallow soil column, coarse vertical discretization, warm bias in soil temperature and lack of phase-dependent thermal conductivity in soil. Furthermore, ERA5L shows a remarkable underestimation of total permafrost area (Table3, Figure 5) when compared with previous estimates. An explanation for this is that the large ALT (i.e. > 1.89 m) that frequently develop in 200 mid-latitude mountains (e.g., Zhao et al., 2010;Cao et al., 2017) cannot be represented by the shallow soil column of ERA5L.
While this could contribute to an underestimation of permafrost area on the QTP, where observed ALT is generally large, we observe a simultaneous cold bias in soil temperature which counteracts the first effect. Our results indicate that the cold bias of ERA5L in mid-low latitudes is highly aligned with the MAAT bias ( Figure 1). This is also suggested by the linear model (Eq. 1). On the other hand, ERA5L underestimates permafrost area in Canada and Alaska despite the observed ALT there being 205 mostly low. This is because the ERA5L soil temperature in western Canada and Alaska is too warm with a wBIAS of about +1.5°C.
Loss of permafrost is an expected consequence of a warming atmosphere. While the loss of near-surface permafrost area derived from ERA5L is similar to previous land-surface model simulations (Lawrence et al., 2008;Slater and Lawrence, 2013), the absolute numbers and the rate of loss have little value for further interpretation because the permafrost area has a 210 pronounced bias to begin with and its temporal dynamics are known to be badly represented with a shallow soil column and are likely affected by an inadequate representation of snow. Furthermore, because permafrost extent is a variable that cannot be observed, we fundamentally lack possibilities for proper validation (Gruber, 2012).

Model evaluation with sparse data
Looking exclusively at summary statistics from 242 sites in 209 grid cells would misleadingly show that ERA5L has a relatively 215 good ability to represent the thermal state of permafrost. For example, consider that ERA5L outperformed the TTOP map in all evaluation metrics (Table 2). However, its simulated permafrost area is visibly low when plotted on a map ( Figure 5).
These contradictory findings can be reconciled because of the warm bias at high latitudes and cold bias in mid-latitudes which cancel each other out based on the available observations ( Figure 3). Clearly, an improvement in summary statistics alone is not a sufficient criterion of superior model performance. Notably, the International Permafrost Association action group 220 "Specification of a Permafrost Reference Product in Succession of the IPA Map" reported in 2016 that, in order to make progress, we needed the capability to measure whether a new map or model output was of superior quality compared with an old one. For this, they recommended that the permafrost community develop and provide the necessary data, methods, and standards (Gruber, 2016).

Scale effects 225
Even for a small area within a single grid cell of Earth-system models or reanalyses (10-100 km), evaluation with point observations remains difficult. This is demonstrated by our simulation example at LdG; within an area of about 20 km × 30 km, MAGT and SO can vary by almost 7°C based on plot sizes on the order of 15 m × 15 m (Gruber et al., 2018). This is important for two reasons. First, the results from statistical evaluations of coarse-scale products such as ERA5L depend significantly on the local selection of observation sites. This issue is known as the spatial effect when the lack of spatially-230 distributed measurements consistent with the size of model grid cells (i.e. 0.1°in ERA5L) is a potential source of error for model evaluation (Gupta et al., 2006;Gubler et al., 2011). Second, ERA5L ground temperatures can only represent at best a small fraction of the area within each individual grid cell. Consequently, their value as part of a permafrost climate services system for informing local decision making (e.g. for adaptation), is limited.

Snow densification and heat transfer 235
The seasonal ERA5L soil temperature deviance ( Figure 2A) and linear model (Eq.1) show a remarkable bias toward high soil temperature in winter that is correlated with snow height. While we do not imply that the GEOtop simulations are correct or accurately represent metamorphism in Arctic snow (see Domine et al., 2019), they do demonstrate that simulations with snow cover of similar mass but different density are able to match ground-temperature observations far better than ERA5L.
Since snow thermal conductivity can be described as an exponential function of its density (Eq. B12), the low-biased snow 240 density of HTESSEL would contribute to a much lower snow thermal conductivity. With the same SWE, a low bias in snow density implies a high bias in snow depth. In this context, the temperature gradient, and hence the heat flux though the snow pack, are both reduced. Using the mean snow density in Table 4 as an example, a snow density of 75% would reduce ground heat-loss though the winter to about 44%. Even though this represents only one local case study at LdG, it sheds light on a possible cause of the ERA5L soil temperature bias in cold regions more broadly. Interestingly, HTESSEL and GEOtop 245 both use the same exponential formulation of snow thermal metamorphism proposed by Anderson (1976) but with different parameters. HTESSEL uses a value of 460 (m 3 kg −1 ) for c ξ , a parameter controlling change in snow density due to thermal metamorphism (Eq. B5) (Dutra et al., 2010). This value is 10 4 times greater than the value for c ξ in GEOtop (Endrizzi et al., 2014) and Anderson (1976). Consequently, for snow densities greater than 150 kg m −3 , the change rate (s −1 ) related to thermal metamorphism remains near zero in HTESSEL. While this may explain, at least in part, the bias in ERA5L snow density and 250 soil temperature, it is unknown whether the excessively high value for HTESSEL is merely an error in the publication cited or whether it reflects the value in the code. An additional contribution of GEOtop to higher snow densities in tundra environments may be the effect of blowing snow (cf, Pomeroy et al., 1993).

Implications
While global reanalyses provide urgently needed meteorological drivers for permafrost simulation, their soil data is not well 255 suited for directly informing permafrost research or local adaptation decisions. As such, simulations using permafrost-specific land-surface models driven by reanalyses (Cao et al., 2019a;Fiddes et al., 2015) will likely be increasingly important in the provision of permafrost climate services. Making future soil-temperature products like ERA5L directly usable will require significant permafrost-specific alterations in model design, especially with respect to snow cover and the total depth of the ground representation for the land-surface models that are used. If indeed the value of the parameter c ξ in the snow metamorphism of 260 HTESSL is in error, then this would be an easy improvement.

Conclusion
Our results support five conclusions.
1 ERA5L soil temperature has a warm-bias at high-latitude and a cold bias in mid-to low-latitude, high-elevation areas.
The soil-temperature bias in permafrost regions correlates with bias in air temperature and with maximum snow height.

265
Seasonally, soil temperatures in winter are more strongly warm biased than in other seasons. With more prevalent snow and ice, ERA5L soil temperatures match observations less well in permafrost-affected regions than in non-permafrost conditions. 4 ERA5L snow density is hypothesized to have a low bias, at least in high-latitude areas, explaining part of the warm bias in soil temperature.
where T obs is observed soil temperature and T mod is the temperature from ERA5-Land soil temperature, GEOtop, or literature.
Appendix B: Snow scheme of HTESSEL

B1 Snow densification
Snow density ρ s (km m −3 ) is constrained to be between 50-450 kg −3 . The compaction of snow density, or change rate (s −1 ), is parametrized as where the first term represents overburden, second term is thermal metamorphism (Anderson, 1976;Boone and Etchevers, 2001), and the last term is the influences of snow liquid water (L s , kg m −2 ) following Lynch-Stieglitz (1994). W s (Pa) is the pressure of overlying snow mass or snow water equivalent (SWE, m) , and η (Pa s −1 ) is the viscosity coefficient of snow.
where g is the acceleration of gravity of 9.807 (m s −2 ). Snow viscosity is described as a function of snow temperature (T s , K) and density following Anderson (1976) where η 0 = 3.7 × 10 7 (Pa s), a η = 0.081 (K −1 ), b η = 0.018 (m kg −3 ). T D (K) is the depression temperature, The change rate of ρ s related to thermal metamorphism is parameterized as where the a ξ , b ξ , and c ξ are constant values of 2.8×10 −6 (s −1 ), 0.042 (-) and 460 (m 3 kg −1 ) derived or modified from Anderson (1976) and Jordan et al. (1999). ∆β s (kg m −3 ) is given as where ρ ξ (kg m −3 ) is equal to 150 kg m −3 . L s is diagnosed from snow temperature, SWE, and snow density, where f (T s ) is the snow temperature function and L c s is the snow liquid water capacity (kg m −2 ).
where T f is 273.16 (K), L c s is parameterized as a function of SWE and β s , where β l s is 200 (kg m −3 ).

B2 Snow cover fraction
Snow cover fraction (SCF) can be given as where SD cr , the minimum snow depth that ensures complete coverage of the grid box, is set as 0.1 m.
B3 Snow thermal conductivity 330 By following Douville et al. (1995), the snow thermal conductivity (λ s ) is treated as a function of snow density, where λ i is ice thermal conductivity of 2.2 (W m −1 K −1 ) and ρ i is ice density of 920 (kg m −3 ). Author contributions. BC carried out this study by analyzing data, performing the simulations, organizing as well as writing the manuscript and was responsible for the compilation and quality control of the observations. SG proposed the initial idea, and contributed to organizing 335 as well as writing the manuscript. XL and DHZ contributed to the writing the paper.
Competing interests. The authors declare that they have no competing interests.
Disclaimer. The authors declare that they have no conflict of interest.
the NPS datasets, thank Xiaoqing Peng for providing the active layer thickness datasets. ERA5-Land reanalysis data and ESA CCI LC map is provided by ECMWF. This study is funded by Strategic Priority Research Program of Chinese Academy of Sciences (grant no.