PIC: Comprehensive R package for permafrost indices computing with daily weather observations and atmospheric forcing over the Qinghai–Tibet Plateau

An R package permafrost indices computing (PIC) was developed, which integrates meteorological observations, remote sensing data, and field measurements to compute the factors or indices of permafrost and seasonal frozen soil. At present, 16 temperature/depth-related indices are integrated into the R package PIC to estimate the possible change trends of frozen soil in the Qinghai–Tibet Plateau (QTP). These indices include the mean annual air temperature, mean annual ground 15 surface temperature, mean annual ground temperature, seasonal thawing/freezing n factor (nt/nf), thawing/freezing degreedays of air and ground surface (DDTa/DDTs/DDFa/DDFs), temperature at the top of the permafrost, active layer thickness, and maximum seasonal freeze depth. The PIC package supports two computational modes, namely, the stations and region calculation that enables statistical analysis and intuitive visualization on the time series and spatial simulations. Over 10 statistical methods were adopted to evaluate these indices in stations, and a sequential Mann-Kendall trend test and spatial 20 trend method were adopted. Multiple visual manners display the temporal and spatial variabilities on the stations and region. The data sets of 52 weather stations and a central region of QTP were prepared and simulated to evaluate the temporal–spatial change trends of permafrost with the climate. Simulation results show extensive permafrost degradation in QTP, and the temporal–spatial trends of the permafrost conditions in QTP were consistent with those of previous studies. The PIC package Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-15 Manuscript under review for journal Geosci. Model Dev. Discussion started: 6 March 2018 c © Author(s) 2018. CC BY 4.0 License.

permafrost table. Changes in permafrost can affect water and heat exchanges, carbon budgets, and natural hazards with the climate change. Permafrost occurs mostly in high latitudes and altitudes with long, cold winters and thin winter snow cover (e.g., the Arctic, Antarctica, Alaska, Northern Russia, Northern Canada, Northern Mongolia, and the Qinghai-Tibet Plateau (QTP)) (Riseborough et al., 2008;Yi et al., 2014b;Zhang et al., 2008a). Over half of the QTP land is underlain by permafrost (Ran et al., 2012). The temperature in QTP has increased by 0.25 °C per decade over the past 50 years (Li et al., 2010;Liu et al., 2006;Shen et al., 2015;Yao et al., 2007). Climate-induced warming of the near-surface atmospheric layer and a corresponding increase in ground temperatures will lead to substantial changes in the water and energy balance of regions underlain by permafrost (Hilbich et al., 2008). Such an increase can warm the ground through energy exchange at the surface and result in significant permafrost degradation. The distribution and changes of permafrost with climate is necessary for infrastructure development, ecological and environmental assessments, and climate system modeling (Luo et al., 2017;Luo et 15 al., 2012;Zhang et al., 2014).
Given the possibility of future climate warming, an evaluation of the magnitude of changes in the ground thermal regime has become desirable to assess the possible eco-environmental responses and their impact on the infrastructure in QTP. Permafrost modeling maximizes quantitative methods, such as analytical, numerical, or empirical methods, to predict the thermal condition of the ground in environments where permafrost may be present (Harris et al., 2009;Lewkowicz and Bonnaventure, 20 2008;Riseborough, 2011;Riseborough et al., 2008;Yi et al., 2014a;Zhang et al., 2008b). At present, dozens of different factors or indices are used to evaluate the characteristics and dynamics of permafrost presence or absence (Riseborough, 2011;Riseborough et al., 2008), including freezing/thawing index, mean annual air temperature (MAAT), mean annual ground Geosci. Model Dev. Discuss., https://doi.org /10.5194/gmd-2018-15 Manuscript under review for journal Geosci. Model Dev. Discussion started: 6 March 2018 c Author(s) 2018. CC BY 4.0 License. temperature (MAGT), mean annual ground surface temperature (MAGST), temperature at the top of permafrost (TTOP), and active layer thickness (ALT), among others. Thereafter, the type and distribution of frozen soil can be classified in a variety of manners depending on the index sizes and magnitude. For example, frozen soil can be divided into highly stable, stable, substable, transitional, unstable, and extremely unstable permafrost, as well as seasonal frozen soil that depends on the size of MAGT (Chen et al., 2012;Ran et al., 2012). These indices can be used to evaluate and predict the temporal and spatial variation 5 in the thermal response of permafrost to the changing climate conditions and properties of Earth's surface and subsurface properties in one, two, or three dimensions (Juliussen and Humlum, 2007;Nelson et al., 1997;Riseborough et al., 2008;Zhang et al., 2005). Accordingly, successfully summarizing and categorizing a variety of frozen soil indices require permafrost modeling that concerns analytical, numerical, and empirical methodologies to compute the past and present condition of permafrost. The Stefan solution (Nelson et al., 1997), Kudryavtsev's approach (Kudryavtsev et al., 1977), TTOP 10 model (Smith and Riseborough, 1996), and Geophysical Institute Permafrost Lab model (Romanovsky and Osterkamp, 1997;Sazonova and Romanovsky, 2003) are several important developments for permafrost modeling in recent years. Permafrost is a subsurface feature that is difficult to directly observe and map. These methods integrate the effects of air and ground temperatures, topography, vegetation, and soil properties to map permafrost spatially and explicitly (Gisnå s et al., 2013;Jafarov et al., 2012;Zhang et al., 2014). Weather observation data, including air and soil temperatures with different depths, are the 15 main inputs for single-point simulation. These indices consist mainly of temperature-related and depth-related indices. The temperature-related indices depict the status of air or land surface temperature in frozen soil environments, whereas the depthrelated indices reveal the status of the active layer. Preparing atmospheric forcing data sets, snow depth and density, vegetation types, and soil classes are generally required for multi-dimensional simulation, which came from multi-source data fusion, particularly remote sensing and ground observation data.

20
The current lack of open source software on permafrost modeling over QTP is a problem. Although many scientists in China have field data and models on hand, the integration of data and models into a new open source model can facilitate the deepening of the discussion and unfolding of permafrost research on QTP. Given the current condition of permafrost modeling in QTP, a comprehensive R package permafrost indices computing (PIC) was developed to compute permafrost and seasonal frozen soil indices. The goal is to determine the solutions to maintain or build the engineers in a manner that provides guidance for the future of highway and high-speed railway design and construction in QTP, as well as further understand the effects of climate change on the permafrost dynamics over QTP. Therefore, the proposed software integrates meteorological observations, remote sensing data, field measurements, and model simulations.
The remainder of this paper is organized as follows. Section 2 describes the prepared data sets, methodology of permafrost modeling, and statistical methods for stations and region. Section 3 presents a detailed description of the functions provided 5 by PIC and the workflow. Section 4 demonstrates the application of the proposed software for the stations and region. Section 5 discusses several benefits and limitations of PIC. Lastly, Section 6 presents the conclusions.

Data and parameters
Daily weather observations. Meteorological data were obtained from the China Meteorological Administration (CMA, 10 http://www.cma.gov.cn/), particularly from permanent meteorological stations across QTP. A total of 52 weather stations with daily meteorological records (i.e., from 1951 to 2010) were selected, including the daily mean, maximum (max) and minimum (min) air temperatures, wind speed, observed and corrected precipitation, evaporation, air humidity, atmospheric pressure, sunshine duration, daily mean, max and min ground surface temperatures, and soil temperature with different depths (i.e., 5, 10, 15, 20, 40, 50, 80, 160, and 320 cm). These data have been corrected under specification for surface meteorological 15 observation and quality control of CMA.
Atmospheric forcing data set. The QT Engineering Corridor (QTEC), which is located at the center of QTP, was selected in preparing the atmospheric forcing data. Global Land Data Assimilation System (GLDAS, https://ldas.gsfc.nasa.gov) and the weather station data of the surrounding QTEC were merged to produce a new data set for 1980 to 2010 with a daily 0.1° temporal-spatial resolution (Rui and Beaudoing, 2011).

20
Parameters. The parameters for the ground conditions were prepared based on vegetation and soil classification (Bicheron et al., 2008;Nachtergaele et al., 2009), field observations, and topographic maps. The parameter data have two data sets: one for weather stations and another for the QTEC region. Table 1 and Figure 1 show the detailed information of the data and parameters.

Permafrost modeling
The PIC package enables the calculation of MAAT, MAGST, MAGT, seasonal thawing/freezing n factor (nt/nf), thawing/freezing degree-days of air and ground surface (DDTa/DDTs/DDFa/DDFs), TTOP, ALT, and maximum seasonal freeze depth (FD). These permafrost and seasonal frozen soil indices that employ the following function were illustrated.

5
As is the annual temperature amplitude at the ground surface, where Tmax and Tmin are the annual max and min temperatures, respectively, at the ground surface. As can be calculated as follows: L is the volumetric latent heat of fusion, ρ is the dry density of soil, and W is the water content of the soil in percentages.
10 DDTa and DDTs are the thawing degree-days of the air and ground surface (Celsius degree-days), respectively. DDFa and DDFs are the freezing degree-days of the air and ground surface (Celsius degree-days), respectively, where Ta and Ts are the air and ground temperatures, respectively, and n is the number of days in a year (Juliussen and Humlum, 2007).
P is assigned a value of 365 days. MAAT and MAGST can be computed as follows:

20
MAGT is the soil temperature in . Tz,t is the ground temperature at any time t and depth z below a ground surface. MAGT can be computed (Juliussen and Humlum, 2007;Riseborough et al., 2008) as follows: nt and nf can be computed (Riseborough et al., 2008) as follows: Two methods serve the same purpose when computing TTOP and ALT. The subscripts S and K stand for the Smith and
The maximum thawing depth or ALT uses the Stefan and Kudryavtsev functions (Kudryavtsev et al., 1977;Riseborough et al., 10 2008), where L is the latent heat of fusion of ice (3.34 × 10 5 J/kg).

15
Freeze_depths is the maximum seasonal freezing depth that uses the Stefan function, which can be computed as follows: Geosci

Statistical methods
Statistical analysis can facilitate the evaluation of the change trend and overall performance of the model simulation. In particular, each statistic has strengths and weaknesses; thus, we adopted over 10 statistical methods to evaluate these indices in station computing for time series data. The quantitative statistics include the slope, y-intercept, Pearson's correlation coefficient (R), coefficient of determination (R 2 ), root mean square error (RMSE), standard deviation (SD), ratio of scatter 5 (RS), normalized RMSE (NRMSE), Nash-Sutchliffe efficiency (NSE), RMSE-observations standard deviation ratio (RSR), percent bias (PBIAS), normalized average error (NAE), variance ratio (VR), and index of agreement (D) (Jafarov et al., 2012;Legates and McCabe, 1999). The sequential Mann-Kendall (MK) trend test was used to statistically assess whether there was a shift in trends of the climate factors and permafrost indices (Fraile, 1993). The original MK trend test can be calculated as follows: Two sequential series ui value can be calculated as follows: The two series for MK trend test, a progressive one and a backward one, were sets up. If they cross each other and diverge 15 beyond a specific threshold value and exceeding the confidence level of 95%, then there is a statistically significant trend shift point.
The spatial trend can also be calculated through the function below. The index represents one permafrost index, n represents the sequential years, and indexi is the index values in year i. Taking ALT as an example, a positive trend value means that ALT was increasing, thereby indicating that permafrost degradation has intensified; a negative value means that ALT was decreasing, 20 thereby indicating that permafrost degradation has a certain inhibition; and a zero trend suggests a lack of change (Chen et al., 2014;Stow et al., 2003). with over 38 functions based on the user's specific requirements (see Figure 2). The following packages are required in setting up the PIC (type library(PIC)): ggplot2 (Wickham et al., 2009), ggmap (Kahle and Wickham, 2013), RNetCDF (Michna and Woods, 2013); (Zambrano-Bigiarini and Rojas, 2013), and animation (Xie, 2013). These packages are automatically added to the users' R library during installation. A data set that contains the daily weather observations, parameters, and information (i.e., from 1951 to 2010) of 52 weather stations in QTP were bundled into this package. However, the region data with the 10 NetCDF format was placed in the GitHub repository. The data set variables excluded in the calculation can be used as reference or provide support to further develop PIC. These variables include wind speed, precipitation, evaporation, humidity, and soil temperature at different depths. PIC was primarily designed to compute indices of permafrost and seasonal frozen soil from observations and forcing data. Therefore, the current stable version of the program (v 1.0) includes functionalities that cover temperature-related indices (i.e., MAAT, MAGST, and TTOP) and depth-related indices (i.e., ALT and FD) that are commonly 15 used in permafrost research.
The PIC package supports two computational modes: the station and region calculations that enable statistical analysis and visual displays on the time series and spatial simulations. The regional calculation adopts GIS approaches to compute each spatial grid. PIC was initially developed to address an immediate need for a reliable and easy-to-use program to estimate the temporal-spatial changes in frozen soil in QTP. Thus, the workflow comprises deliberately simplified steps involved 20 throughout the entire process. Once PIC is installed, the workflow of the weather observations is considerably straightforward: (1) an index of a weather station for one year or multiple years is calculated, (2) an index of 52 weather stations from 1951 to 2010 is calculated, and (3) an index of all stations or permafrost stations from 1951 to 2010 is drawn through curve and spatial visualization.
Step (1) is an optional step. The workflow of the forcing data has only two steps: (1) a total of 4 indices from  (2) the spatial statistic and visualization of these 4 indices are drawn. Table 2 describes most of these functions.

Examples
Several examples of the PIC use and application were presented here. This section highlights several significant features of the package in terms of specific functions, including station and region calculation, statistics, and visualization. However, PIC VarName option should add "_air" or "_ground" at the end of the Freezing_index and Thawing_index. However, the abbreviation can also be utilized as the option input. The "Thawing_index_air" and "ta" are the same.
Com_Indices_QTP ( A spatial trend can also be computed using the Spatial_Stat function after the regional calculation. The function simultaneously saves the spatial trend of the five indices into the NetCDF file. In addition, the function draws the animation of the spatial trend (see Section 4.4).

Visualization
Station visualization can be produced by Plot_TTOP_ALT and Plot_3M. The Plot_TTOP_ALT function plots two TTOP or two ALT indices in a figure for all stations or stations with permafrost. VarName has the "TTOP" and "ALT" options, whereas SID has the "permafrost" and "all" options. The Plot_3M function draws the MAAT, MAGST, and MAGT indices. The two functions only plot these stations where permafrost exists when SID = "permafrost."

Map_Pic (VarName = "TTOP_S")
Map_Pic (VarName = "TTOP_K") The input and output of the regional calculation can be drawn using the Netcdf_Multiplot function (see Figure 6); the Netcdf_Animation function uses animation to display these values. The spatial trend can also be drawn in the Spatial_Stat apart from calculating the spatial statistics. This function draws all four indices when "VarName" has no input (see Figure 7). The PIC package computes and maps the temporal dynamics and spatial distribution of permafrost in the stations and region.
The regional modeling underwent more challenges than the stations' input data and parameters. The station calculation can estimate the long-term temporal trend of the permafrost dynamics, whereas regional calculation can estimate the temporal-15 spatial trend. Climate change indicates a pronounced warming and permafrost degradation in QTP (Chen et al., 2013;Cheng and Wu, 2007;. The simulation results show widespread permafrost degradation in QTP, and the temporal-spatial trends of the permafrost conditions in QTP were consistent with those in previous studies . In addition, the simulated TTOP and ALT that uses the Stefan and Smith functions have higher TTOP and ALT than the Kudryavtsev function. Although the overall trend of TTOP and AIT are coincidental, two different computational methods can be combined to simulate their variations. Furthermore, 16 indices can be collectively employed for a comprehensive analysis. Moreover, the station and region modeling can be integrated to evaluate the temporal-spatial 5 evolution of permafrost in the QTP. In particular, the station modeling can be applied to validate the simulated results of the region. Moreover, regional calculation can extend from QTEC to the entire QTP.
The "for" loop is discarded, whereas the "apply" functions are used extensively to significantly lower the computation time.
The current regional calculation only takes approximately 11 s. Apart from the Kudryavtsev model that requires considerable computation time (i.e., approximately 5 min), the station calculation also exhibited an improved efficiency. Therefore, PIC 10 can be considered an efficient R package.

Advantages
Previous studies in the QTP (1) used one or two indices, such as MAAT and MAGST, to evaluate the permafrost changes ; (2) presented a static permafrost distribution that constructed a regression analysis method through the relationship between MAGT and elevation, latitude, and slope-aspects (Nan, 2005;Yin et al., 2017); and (3)  weather stations of QTP were introduced, which can approximately resolve the spatial change directions of the permafrost area. The QTEC region is an example of spatial modeling, which classifies land cover and topographic features to determine the input spatial parameters. Spatial modeling also uses the GLDAS satellite data to provide detailed information of the active layer and permafrost. The static/dynamic maps and statistical values of these indices can facilitate the understanding of the current condition of the near-surface permafrost and identify stations and ranges at considerably high risk of permafrost 5 thawing with the changing climate and human activities. Permafrost thawing causes significant changes in the environment and characteristics of frozen-soil engineering (Larsen et al., 2008;Niu et al., 2016). A comprehensive assessment of permafrost can provide guidance regarding the future of highways and high-speed railway systems in QTP.

Limitations and uncertainties
PIC was developed with numerous indices, as well as support station and regional simulation. The PIC package can be used 10 to estimate the frozen soil status and possible changes over QTP by calculating permafrost indices. This package has many engineering applications and can be used to assess the impact of climate change on permafrost. Moreover, this package provides observation data and provides the ability of comprehensive analysis through multiple indices. The probability of permafrost occurrence and most likely permafrost conditions are determined by computing the 16 indices. Although PIC quantitatively integrates most of these indices based on previous studies (Jafarov et al., 2012;Nelson et al., 1997;Riseborough 15 et al., 2008;Smith and Riseborough, 2010;Zhang et al., 2005;Zhang et al., 2014), it still has several limitations.
First, the regional calculation is one-dimensional and assumes that each grid cell is uniform without the water-heat exchange.
Second, soil moisture changes at different depths affect the thermal conductivity and thermal capacity of the soil (Shanley and Chalmers, 1999;Yi et al., 2007); thus, the soil input parameters should be dynamically changed. Lastly, climate forcing has several uncertainties , including input air and ground temperatures (i.e., the quality of the ground 20 temperature in QTP is currently unreliable); thus, the regional calculation supports fewer indices than the station calculation.
These deficiencies can be significant for the permafrost dynamics with environmental evolution.

Conclusions
An R package PIC that computes the temperature/depth-related permafrost indices with daily weather observations and climate forcing has been developed. This package is open source software and can be easily used with input data and parameters, and users can customize their own data and parameters. A total of 16 permafrost indices for stations and region are currently integrated into the R package PIC to estimate the status of the active layer and permafrost in QTP. The current functionalities 5 also include time-series statistics, spatial statistics, and visualization. Multiple visual manners display the temporal and spatial variability on the stations and region. The package produces high-quality graphics that illustrate the status of frozen soil and may be used for subsequent publication in scientific journals and reports. The data sets of the 52 weather stations and a central region of QTP were prepared and simulated to evaluate the temporal-spatial change trends of permafrost with the climate. The simulated PIC results generally indicate that the temporal-spatial trends of permafrost conditions essentially agree with 10 previously published studies. The PIC package has many engineering applications and can be used to assess the impact of climate change on permafrost. Additional features may be implemented in future releases of PIC to broaden its application range. In the future, the observation data of the active layer will be integrated into the PIC data set, and the output will be compared with the observation data. The PIC package can also be used to predict the future state of permafrost by utilizing projected climate forcing and scenarios. Additional functions and models will be absorbed into PIC to improve the simulation 15 performance and perform comparative analyses with other functions and models. Parallel computation will be added to the PIC package to improve the computation efficiency. The key impact that PIC is expected to provide on the open community is an increase in consistency within and comparability among studies. Furthermore, we encourage contributions from other scientists and developers, including suggestions and assistance, to modify and improve the proposed PIC package.

20
The PIC code that support the findings of this study is stored within the GitHub repository (https://github.com/iffylaw/PIC).

Data availability
The data is included in the Supplement files or GitHub repository.

Competing interests
The authors declare no competing financial interests.