Prediction of Seasonal Frost Heave Behavior in Unsaturated Soil in Northeastern China Using Interactive Factor Analysis with Split-Plot Experiments and GRNN

: Few studies of frost heave mechanisms have considered multifactor interactions, particularly in unsaturated saline soils typical of northeastern China. We collected soil samples in western Jijin Province and assessed their potential frost heave behavior with reference to four controllable factors: soluble salt content (CSS), compactness (C), temperature (T), and water content (WC) using a two-level split-plot experiment. The resulting frost heave ratio was between − 0.6% and 2.1%. Analysis of variance showed that water content, compactness, and temperature had signiﬁcant e ﬀ ects on frost heave behavior, with water content having the strongest correlation (factor coe ﬃ cient of 0.82), while content of soluble salt (CSS) had no signiﬁcant e ﬀ ect. The interaction factors (products of single factors) CSS × WC and C × WC had signiﬁcant e ﬀ ects on frost heave behavior. A correlation analysis using these interaction factors with experimental data drawn from previous research showed results consistent with the improved frost heave experiment as the signiﬁcant e ﬀ ects of single factors on frost heave behavior ranked from WC > C > T and the interaction factors CSS × WC and C × WC gain had signiﬁcant e ﬀ ects. We then established two generalized regression neural network (GRNN) models based on the single and interaction factors in order to predict frost heave behavior, showing that adding the latter to the input dataset improved the model accuracy. Thus, future research on predicting frost heave behavior in unsaturated saline soils should consider multiple interacting factor for greater accuracy. Contributions: Conceptualization, H.W. and J.B.; Data curation, H.W.; Formal analysis, H.W.; Funding acquisition, J.B. and Y.W.; Methodology, H.W. and J.W.; Project administration, H.W. and X.S.; Resources, X.S.; Validation, H.W. and J.W.; Visualization, H.W. and Z.J.; Writing—original draft, H.W.; Writing—review & editing, H.W. and Y.W.


Introduction
Frozen soil can be divided into permafrost, seasonal frozen soil, and short-term frozen soil [1]. In cold regions, freeze damage caused by frozen soils seriously threatens ecological functions and engineering safety. Unsaturated soil is widely distributed in cold regions; when the temperature drops below the freezing point, the water in such soil begins to freeze [2]. During this process, soil water, vapor, and heat move from warm to cold locations, driven by a complicated process involving the temperature gradient and soil water potential. At the same time, discontinuous ice lens growth in the unsaturated soil causes frost-heaving [3][4][5], leading to soil body expansion and uneven uplift of the ground surface [6]. This process can cause serious damage to above-ground and underground infrastructure, resulting in serious economic losses [1, 7,8]. For example, uplift and cracking of road surfaces or changes in railway grades can threaten safe transportation [9] or disrupt buried pipes

Study Area
Northeast China experiences large-scale seasonal freezing and is widely characterized by saline soils [28]. The study area in western Jilin Province has a climate that is dry and windy in spring, hot in summer, and cold in winter. Evaporation tends to be far greater than precipitation, resulting in salt accumulation at the soil's surface. Daily mean temperatures range from 38 • C to −36.5 • C with an annual average temperature of 4.6 • C. The frozen soil period lasts from November to March, with freeze depth reaching 80-100 cm. As the saline soil begins to form a frozen soil layer, the soil water moves from the unfrozen layer to the frozen layer under the influence of total water potential, carrying solutes to the surface and exacerbating soil salinization [29]. Increased water content at the soil's surface promotes more intense frost heave with associated risks to infrastructure.

Sample Collection
The saline soil used in the experiment was collected from Da'an city in western Jilin Province during May 2009. Soil samples from different depths (20,30,40,50,70, 100 cm) were collected at a single sampling point using the cutting ring method (Φ100). As the disturbed soil samples would inevitably lose water during transportation to the laboratory, smalls amount (~20-30 g) of original soil (3 samples per layer) were collected in aluminum cases, weighed on-site, and kept sealed; these were used to calculate water content while other original samples were used to test soil density. Disturbed soil samples were also collected at each layer, each weighing at least 1000 g; these were used to test granulometric composition, content of soluble salt (CSS), and pH, as well as conducting the improved frost-heaving experiments. Examples of original and disturbed soil samples are shown in Figure 2.

Soil Testing
All soil samples were brought back to the laboratory on the same day. Disturbed soil samples were naturally dried, ground, sieved (2 mm), and stored. Granulometric composition, water content, content of soluble salt (CSS), pH, and density were tested at the Key Laboratory of Jilin University. A laser particle size analyzer (Mastersizer 2000) was used to analyze the granulometric composition in terms of clay (<0.002 mm), silt (0.02-0.002 mm), and sand (2-0.02 mm). The remaining indicators were tested following the methods introduced by Bao [30]. All test results are given in Table 1. Water content (density) values are the mean of the 3 original samples for each layer. Granulometric composition, CSS, and pH values reflect the disturbed samples.
Soil samples at different depths clearly had similar granulometric compositions. The fine-grained soil (clay and silt) content values were all high, defining all soil samples as silty clay. Water content was highest at 40, 50, and 70 cm, while the shallowest sample had the lowest water content. The CSS values decreased consistently with depth. As soil with high water content and CSS is most likely to cause serious frost heave activity, we chose to use disturbed soil samples at 30 cm depth for further experimentation.

Improved Frost Heave Experiments
Unlike past single-factor studies, our improved frost heave experiments considered interactions between four factors: soluble salt content (CSS), compactness (C), temperature (T), and water content (WC) ( Table 2), where the interaction factors are the products of two or three single factors. We used disturbed soil samples at 30 cm depth and determined factor levels with reference to soil properties and experimental data from previous research [22]. CSS was controlled by adding NaHCO 3 to soil samples because the study area was composed of carbonate saline soil. Water content was controlled by adding water [22]. A total of 16 experiments were conducted, divided by temperature conditions (either −25 • C or −4 • C) into two groups that each lasted 12 hours at a fixed temperature (Table 3).  Unsaturated soil was placed in an experimental tube with an inner diameter of 5 cm and a height of 10 cm; the soil column height was 8 cm. The tubes were made of thermal insulation material and their bottoms were sealed with the same material to ensure that the soil gradually froze from top to bottom during the experiment, similar to actual natural conditions. A dial indicator was placed above each tube to measure the longitudinal volume change of the soil column ( Figure 3). The prepared samples were immediately placed in a temperature control box. We defined the frost heave effect η by the following ratio.
where ∆H (cm) is the height increment of the soil column after freezing and H 0 (cm) is the original height of the soil.

Frost heave Experiments
Bao [22] collected soil samples from the study area and carried out 810 frost heave experiments analyzing the individual influence of same four factors as in our study (Table 4). However, as that study used the OFAT approach, the interactions between these factors could not be determined. Table 4. Influence levels of four factors in frost heave experiments [22].

Factors Levels
Compactness (

Data Analysis
We conducted all data analysis, including restricted maximum likelihood (REML), F-test, data diagnostics, and graphing, using Design-Expert 10 (DX10) [31] software. REML is a method for estimating the parameters of a linear regression model. The F-value associated with the regression model is estimated based on the techniques proposed by Kenward and Roger (the KR method) [32]. This method uses the components of variance as estimated by REML, the model being fit, and the design structure to approximate the correct distribution for testing the terms and model subsets. Data diagnostics and graphs can help determine the prediction accuracy of the model (whether predicted values were close to experimental values) [33]. R 2 and adjust-R 2 (Adj-R 2 ) were used for quantitative assessment of the regression model's prediction accuracy. The former reflects the accuracy of the regression model by showing the degree to which the input variables explain the variation of the output variable. The latter is similar, but also accounts for the number of variables. For the experimental data of Bao [22], we added interaction factors and used correlation analysis between the response (η) and influencing factors.

Prediction
We used GRNN to predict frost heave behavior using the four influencing factors as input data and the frost heave ratio as the output target. All 810 experimental data points from Bao [22] were randomly divided into 3 groups: a learning dataset (680 data points) for GRNN, a training dataset (100 data points), and a test dataset (30 data points). In order to analyze the influence of interaction factors on the frost heave prediction accuracy, we prepared two types of input datasets: single factors only (Single-GRNN) and single factors with interaction factors CSS×WC and C×WC (Interact-GRNN), using R 2 and root mean square error (RMSE) to evaluate model accuracy. GRNN is a radial basis function neural network proposed by Specht [34] that has a simple training procedure and unique adjustment parameters. The GRNN network structure is divided into four layers (input, pattern, summation, and output; Figure 4). The input vector is X = [x 1 , x 2 , ..., x m ] T and the corresponding output is Y. The dimension of the input vector X is m and the number of samples is n. This model can automatically adjust the network structure according to the training samples; the number of nodes need not be determined by the user. The nodes of input layer are merely distribution nodes, the number of input nodes is equal to the dimension of the input vector, m. Arrows in Figure 3 indicate data input, transfer, and output. The input layer provides all X to the second layer, the pattern layer. The number of nodes on the pattern layer is equal to the number of samples, n. The activation function is given in Equation (2).
where p i is the activation function of i-th node, X is the input vector, X i is the learning sample of i-th node, δ is the smoothing factor, and upscript T means matrix transposition. The pattern layer outputs are passed to the summation layer. Two ways of summation are performed in summation layer. One is f (X)K which is the summation of the outputs of the pattern layer (Equation (3)), the other isŶ f (X)K, which is the summation of the outputs of the pattern layer weighted by the observation output variable Y i (Equation (4)).
where K is a constant that does not need to be computed. The output layer dividesŶ f (X)K by f (X)K to yield the desired estimate ofŶ(X), see Equation (5).

Results of Improved Frost Heave Experiments
The resulting η values ranged from −0.60% to 2.1% ( Table 3). Five of the sixteen response values were negative, indicating that the saline soil shrank rather than swelled under some low-temperature conditions. The coefficient estimates, F-value, and p-value showed that factors compactness, temperature, and water content all had significant effects on η (p-value < 0.05) ( Table 5). The effect of CSS on η was not significant, but the interaction factor CSS × WC had a significant effect on η, indicating that CSS alone had little effect on frost heave behavior but did interact with water content to create an effect. The interaction factor C × WC also had a significant effect on η. Water content had the largest coefficient of interaction (0.82), indicating that the frost heave ratio increased with higher water content; the next-highest factors were CSS × WC (0.21) and C × WC (−0.18). Temperature had the smallest coefficient (−0.078), indicating that frost heave ratio increased only slightly with a decrease in temperature at low-temperature conditions (below 0 • C).
The interactions between CSS and water content, and compactness and water content, are shown in Figure 5, where the two intersecting lines indicate that the effect of one factor depends on the level of the other. The slope of the line can be used to indicate the sensitivity of η to variables. For example, at low CSS (high compactness), η is less sensitive to water content than at high CSS (low compactness). Thus, if we wanted to make the frost heave ratio even more robust to variations in water content, we could reduce CSS (increase compactness). The liner regression model Equation (6) was set up using the factors ranked in descending order of absolute coefficient value (Table 5) to develop the relationship between experimental and predicted results for η ( Figure 6). The R 2 of the model was 0.9951 and the Adj-R 2 was 0.9754, indicating that the predicted values matched the experimental data. The 3D response surface plots for the effects of CSS and compactness on η at different levels of temperature and water content ( Figure 7) show that η had an inverse relationship with temperature, which was most obvious when water content was high (26%) than low (18%). For a given amount of water, lower temperature produced a more complete phase change from liquid to solid such that an increase in ice content can lead to increased η [22]. In contrast, the relationship between η and CSS changed at different water content levels; η increased with increasing CSS when water content was high, but decreased with increasing CSS when water content was low. Compactness and η had a positive relationship when CSS was between~0.5-0.7%, however when CSS was higher, this effect vanished or even turned inverse. Water content had a clearly positive relationship with η; when more soil water is available for freezing, more ice lenses form and cause frost-heaving.  For single factors, the results of the improved frost heave experiments indicated that the frost heave ratio was indeed a function of water content, compactness, and temperature, consistent with conclusions obtained from previous OFAT trials [35]. Sheng et al. [13] analyzed the sensitivity of the frost heave ratio of clay to compactness, temperature, and groundwater level, showing that this increased significantly as groundwater level increased, causing the water content of clay to increase and providing source material for ice lens formation and aggravated frost heave [36,37]. This demonstrates the positive correlation between water content and frost heave ratio. Similarly, the effect of temperature on frost heave in this study was consistent with past research [9,13,17].
Compactness was positively correlated with frost heave in this study. Uncompacted soil has high pore space that frozen ice lens cannot fill, while compacted soil has tightly structured soil particles with low pore space that is more easily filled by frozen ice lenses, leading to obvious frost heaving deformation [38]. This result was consistent with Sun et al. [21] and Bao [22], but other studies have concluded that increasing compactness reduces frost heave [13,39,40]. This contrast shows that the relationship between compactness and frost heave is affected by interactions [23] between compactness and other factors such as water content, as demonstrated in our analysis. The interaction factors CSS × WC and C × WC showed that these dynamics represent the failure of one factor to produce the same effect on the response at different levels of another factor [23]. When optimizing responses (minimize/maximum η) during engineering design and construction, it is necessary to consider such interactions. For example, in soil-filling projects (such as roadbeds, embankments, or dams), using low-CSS soil or mixed soil (adding coal ash or other materials) to fill the foundation, laying an anti-seepage layer to block moisture exchange, and increasing compactness all help to reduce the risk of frost heave [41][42][43].

Frost Heave Experiments
Soil parameter statistics and experimental results are compared in Table 6; 810 experiments were conducted by Bao [22] and complete data are given in Appendix A (Table A1). In order to consider the influence of interaction factors on η, CSS × WC and C × WC were added to the experimental data; the results of correlation analysis using SPSS 22.0 software are shown in Table 7. The three correlation tests had consistent conclusions. The absolute value of the correlation coefficient between C × WC and η was largest (significant at the 0.01 level). CSS × WC was also significantly related to η at the 0.01 level. Of the four single factors, water content had the largest absolute coefficient value (positive), indicating that water content had a significant positive correlation relationship with η. Compactness was also significantly positively correlated with η, while temperature was negatively correlated. The Pearson test showed that the correlation between CSS and η was weakly positive at the 0.05 level, but the Kendall and Spearman test showed that the correlation between CSS and η was not significant. These results were consistent with the 16 improved frost heave experiments in showing that CSS had no significant effect on the frost heave ratio, while the others increased in effect from T < C < WC; the interaction factors CSS × WC and C × WC both had significant effects.
We used experimental data (Table A1) for correlation analysis to determine the relationship between frost heave ratio and different factors (consistent with Bao [22]) and the influence degree of factors on frost heave by correlation coefficient. After the interaction factors were added to the experimental data, the subsequent correlation analysis showed that the interaction factors did have a significant effect on frost heave, supporting our initial results. The conclusions from our 16 improved frost heave experiments were similar to those of the 810 OFAT frost heave experiments from Bao [22]: workload was reduced by 98.02%, demonstrating that the split-plot experimental design was helpful for obtaining reasonable results with reduced cost and workload. More importantly, interaction factors were more properly taken into account in the split-plot frost heave experimental design. Table 8 shows the statistical indexes of the GRNN models. Figure 8 shows the results of the GRNN models' training and testing.  The training and testing R 2 values of the two GRNN models were all close to 1, indicating that the predicted values were very close to the experimental observations and demonstrating that the GRNN method can accurately predict frost heave behavior in the saline soil from the study area. Although the predicted values for Single-GRNN deviated from the 1:1 line, those for Interact-GRNN were all near this line, suggesting that the latter's prediction accuracy was higher. Table 8 supports this conclusion, as the training and testing R 2 values for Single-GRNN were 0.95 and 0.92, respectively, compared to those for Interact-GRNN (0.99 and 0.97, respectively), an increase of 3.62% and 5.14%, respectively. Similarly, the Single-GRNN training and testing RMSE values were 0.2000 and 0.2370, compared to those for Interact-GRNN (0.1059 and 0.1373, respectively), a decrease of 47.02% and 42.08%, respectively. Thus, we conclude that including interaction factors to the input dataset can significantly improve the GRNN results when simulated frost heave behavior of saline soil by reducing the deviation between predicted and observed values. In comparison, Zhang et al. [20] used BPNN and GRNN to obtain the frost heave ratio from data for water content, compactness, temperature, and content of soluble salt. Those results demonstrated that both methods were suitable for predicting the frost heave ratio of saline soil although the GRNN model was more accurate. However, that study did not include interaction factors in the input dataset for either model to improve the outputs. The testing R 2 (RMSE) value of GRNN in Zhang et al. [20] is 0.87 (0.2598), which is smaller (bigger) than R 2 0.97 (RMSE, 0.1373) of Interact-GRNN.

Conclusions
In order to determine interaction factors affecting frost heave and establish a high-accuracy GRNN model to predict the frost heave behavior of saline soil, split-plot frost heave experiments were conducted and a GRNN model including interaction factors was applied, with the following conclusions: Of the single factors, water content had the greatest influence on the frost heave ratio (η) of unsaturated soil with a significant positively correlation. Factor compactness was also significantly positively correlated with η, while temperature was negatively correlated with η and CSS had no significant correlation. While these results are applicable to the study area's sodic saline soil and are consistent with other research [21,22], some studies on frost heave have shown that compactness is negatively correlated with frost heave. This is because the interactions between compactness and water content affect the relationship between compactness and frost heave. The split-plot experiment was able to analyze the influence degree of all single factors while considering the relationships between interaction factors and the response value. This study design can also be used in other multifactorial experiments including interaction factors and significantly reduces the experimental workload while optimizing the experimental arrangements. We only set two levels (low and high) for each factor and used a linear model to represent the relationship between frost heave and the influencing factors, but future research could set more factor levels and use a nonlinear model. Our results showed that interaction factors CSS × WC and C × WC had significant effects on the frost heave behavior of saline soil, as supported by correlation analysis. At low CSS (high compactness), the sensitivity of η was smaller than at high CSS (low compactness). This could be important for engineering design and construction in cold regions: for example, foundations could be covered with low-CSS soil to increase compactness and control frost heave damage. The accuracy of our results demonstrated that the GRNN model was suitable for predicting frost heave ratios in unsaturated saline soil in the study area, and that its performance was improved by the inclusion of interaction factors rather than single factors only.

Acknowledgments:
The authors wish to thank Shuochao Bao for providing complete data of single-factor frost-heaving experiment.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.