Development of Pedo-Transfer Functions for the Saturated Hydraulic Conductivity of Forest Soil in South Korea Considering Forest Stand and Site Characteristics

: The saturated hydraulic conductivity ( K s ) is one of the most important soil properties for many hydrological simulation models. Especially in South Korea, analyzing the K s of the forest soil is essential for understanding the water cycle throughout the country, because forests cover almost two-thirds of the whole country. However, few studies have focused on the forest soil in the temperate climate zone on a nationwide scale. In this study, 1456 forest soil samples were collected throughout South Korea and pedo-transfer functions employed to predict the K s were developed. The non-linearities of the soil and topographic features were considered with the pretreatment of variables, and the variance inﬂation factor was used for treating the multicollinearity problem. The forest stand and site characteristics were also categorized by an ANOVA and post hoc test due to their diversity. As a result, the K s values were di ﬀ erent for various forest stands and site characteristics, which was statistically signiﬁcant. Additionally, the model performance was higher when both soil properties and topographic features were considered. The sensitivity analysis showed that the K s was highly a ﬀ ected by the bulk density, sand fraction, slope, and upper catchment area. Therefore, the topographic features were as important in predicting the K s as the soil properties of the forest soil.


Introduction
The saturated hydraulic conductivity (K s ) is an important factor that represents the basic properties of soil. It can represent the rate of infiltration, so it is essential for understanding the water cycle through soil, such as water recharge, drainage, baseflow, and runoff generation [1]. To date, many hydraulic simulation models, such as TOPMODEL, HYDRUS, and DHSVM, have been developed to simulate the water flow in a catchment [2][3][4]. These models have several input variables for simulation, and the K s is included in the model equations, which means that the K s can directly influence the model outcomes [4]. Therefore, the K s , which is relevant to the infiltration rate, is one of the most important input factors and highly affects the model output, such as runoff [5,6]. Moreover, in order to achieve the Sustainable Development Goals (SDGs) related to ensuring the availability and sustainable management of water, the estimation of soil properties that are closely related to the water yield simulation and forest water management is important [7,8].

Soil Physical Properties and Topographic Features
Forest soil samples were collected through a 100 cc soil core sampler, and a total of 1456 soil samples were collected throughout the forest nationwide. Six properties were analyzed through laboratory experiments to investigate the physical properties of the soil. The analyzed soil properties were the saturated hydraulic conductivity; bulk density; sand, silt, and clay fraction; and organic matter ( Table 1). The distribution of the 1456 soils across the USDA (United States Department of Agriculture) textural classes and the range of saturated hydraulic conductivity ( ) by textural class are shown in Figure 2.

Soil Physical Properties and Topographic Features
Forest soil samples were collected through a 100 cc soil core sampler, and a total of 1456 soil samples were collected throughout the forest nationwide. Six properties were analyzed through laboratory experiments to investigate the physical properties of the soil. The analyzed soil properties were the saturated hydraulic conductivity; bulk density; sand, silt, and clay fraction; and organic matter ( Table 1). The distribution of the 1456 soils across the USDA (United States Department of Agriculture) textural classes and the range of saturated hydraulic conductivity (K s ) by textural class are shown in Figure 2.   The elevation, slope, topographic wetness index (TWI; Equation (2)), upper catchment area (CA), plan curvature, and profile curvature were collected as the topographic features. The elevation and slope indicate the altitude and local slope at the site where the soil was collected, respectively. The topographic wetness index, also known as the compound topographic index (CTI), indicates the steady state wetness index by topographic features. The TWI is calculated by the local upper catchment area and local slope and is derived as follows: where a is the local upper catchment area draining through a certain point per unit contour length and is the local slope [20]. The plan curvature can represent the shape of the horizontal plane, and the profile curvature can represent the shape of the vertical plane. The TWI, plan curvature, and profile curvature were calculated with the GIS spatial analyst tool.
In addition, the characteristics of the forest stand and site from which the soil samples were collected were investigated. The forest stand was divided into three categories: coniferous forest, broadleaf forest, and mixed forest. Bedrock and landform were investigated as site characteristics. Bedrock was largely divided into five types, including igneous, sedimentary, and metamorphic rock, and landform was classified into five types, as defined in Table 2. The elevation, slope, topographic wetness index (TWI; Equation (2)), upper catchment area (CA), plan curvature, and profile curvature were collected as the topographic features. The elevation and slope indicate the altitude and local slope at the site where the soil was collected, respectively. The topographic wetness index, also known as the compound topographic index (CTI), indicates the steady state wetness index by topographic features. The TWI is calculated by the local upper catchment area and local slope and is derived as follows: where a is the local upper catchment area draining through a certain point per unit contour length and tanb is the local slope [20]. The plan curvature can represent the shape of the horizontal plane, and the profile curvature can represent the shape of the vertical plane. The TWI, plan curvature, and profile curvature were calculated with the GIS spatial analyst tool. In addition, the characteristics of the forest stand and site from which the soil samples were collected were investigated. The forest stand was divided into three categories: coniferous forest, broadleaf forest, and mixed forest. Bedrock and landform were investigated as site characteristics. Bedrock was largely divided into five types, including igneous, sedimentary, and metamorphic rock, and landform was classified into five types, as defined in Table 2. We selected PTFs that have frequently been cited in the literature for comparing the K s of forest soils in South Korea to the K s from other international databases by estimating the saturated hydraulic conductivity using our database. Six PTF models were selected, and brief explanations and references are shown in Table 3. Subsequently, the forest soil characteristics of South Korea were compared with other soil conditions by comparing the measured and predicted saturated hydraulic conductivity, which were estimated by six PTFs with the measured soil physical characteristics as input data. Table 3. Six pedo-transfer functions (PTFs; Equations (S1)-(S6)) of previous research used for estimating the saturated hydraulic conductivity and their required soil properties.

K s Unit Required Soil Properties for Estimation Land Use Types
Puckett et al. [12] mm h −1 Clay The lower coastal plain Campbell and Shiozawa [21] mm h −1 Sand, Clay na Jarbro [11] cm day −1 Silt, Clay, ρ b na Dane and Puckett [17] mm h −1 Clay The lower coastal plain Wosten et al. [15] cm day −1 Silt, Clay, Organic matter, ρ b na Julia et al. [14] mm h −1 Sand, Clay, Organic matter na Note: ρ b is the bulk density and na is not available. Soil databases from many previous studies were based on various land use types, so there was no clear information about land use types.
Six PTFs did not focus on the forest soil, and many previous studies were based on various land use types. There was also no clear information about the land use types because previous research was conducted primarily focusing on local regions in which soil surveys were carried out. The PTFs in previous research have been empirical equations. Therefore, the ranges of the soil properties from which each PTF was derived are represented in Table 4.

Preprocessing of Explanatory Variables
In this study, PTFs were developed with multiple linear regression analysis to estimate the logarithmized saturated hydraulic conductivity (ln(K s )). Since the saturated hydraulic conductivity can have a non-linear relationship with soil characteristics and topographic features, logarithmized (ln(x)), exponential (e x ), squared (x 2 ), and untreated (x) treatments were carried out, and linear correlations with ln(K s ) were compared. Subsequently, the most linearly correlated variable type was adopted as the explanatory variable for multiple linear regression. The Pearson correlation coefficient was used to confirm the linear correlation.

Detecting Multicollinearity Using the Variance Inflation Factor (VIF)
In multiple linear regression, when there is a correlation between variables, coefficients in the multiple regression equation can be unreasonable and unreliable values, which is called a multicollinearity problem. Because this happens when a linear correlation occurs between the explanatory variables, it is necessary to remove the highly correlated variable in order to produce a reasonable result [22]. The variance inflation factor (VIF) is one of the best methods for identifying the correlation between variables. Since the VIF is a numerical value after multiple linear regression analysis is performed between explanatory variables, the higher the linear correlation between the explanatory variables, the higher the value of the VIF. If the VIF is higher than 5, it is determined that there is multicollinearity and the variable is removed [23]. The variance inflation factor (VIF) for one explanatory variable (x 1 ) is derived as follows: where R 2 is the coefficient of determination of the multiple regression equation, which is estimated with one variable x 1 as a response variable and other variables as explanatory variables.

Model Assessment
Seventy percent of the total data set was used for developing the PTFs, and the remaining 30%, the hold-out test data set, was used to verify the model. To assess the performance of the PTF model developed in each process, the root mean squared log-transformed error (RMSLE), mean log-transformed error (MLE), and coefficient of determination (R 2 ) were used, and these are defined as follows: where n is the number of observation samples, K s is the measured saturated hydraulic conductivity, K s is the predicted value from the PTFs, and K s is the mean of the measured values. The RMSLE and MLE represent the differences between the measured and predicted values as absolute values, and the lower the performance of the model, the higher the values. The RMSLE is directly affected by the scale factor and reacts more sensitively to outliers than the MLE. R 2 is the dimensionless value that represents the correspondence between the measured data and the predicted data. This has a value between 0 and 1, and variables are more correlated as the R 2 nears 1.

Sensitivity Analysis
Sensitivity analysis can determine the effect of each explanatory variable on the response variable by confirming the amount of change in the response variable as the explanatory variable changes [24]. The coefficient of multiple linear regression (MLR) can represent the sensitivity of the variable because the explanatory and response variables have a linear relationship. Therefore, it is not necessary to conduct sensitivity analysis if MLR is adopted as a model. However, in this study, several MLRs were presented, not just one MLR, to estimate the saturated hydraulic conductivity. Additionally, each explanatory variable had nonlinearity throughout the preprocessing.
Each explanatory variable was varied by total 10 multipliers: 0.1, 0.25, 0.5, 0.8, 0.9, 1.1, 1.2, 2, 4, and 10. Each variable was multiplied while the others were kept constant. Moreover, the multiplied variable was used as the input variable in PTFs, and the output, which was the modified saturated hydraulic conductivity, was normalized (modified ln(K s )/non-treated ln(K s )) to determine the amount of change. A sensitivity run was only implemented once for the PTF model, because there was no variance between the output of each sensitivity run.
Statistical analyses, such as the Pearson correlation analysis, analysis of variation (ANOVA), multiple linear regression, and p-value determination, were conducted using the SPSS software. We confirmed that the result was statistically significant and rejected the null hypothesis when the p-value was less than 0.05.

Saturated Hydraulic Conductivity in Forest Soil
Six PTFs were analyzed to confirm the difference between the K s of the forest soil collected in this study and the K s predicted by the PTFs in previous research. The K s was calculated using 1456 soil physical characteristics as input factors for the PTFs in previous research. Six PTFs did not focus on the forest soil, and many previous studies were based on various land use types. There was also no clear information about land use types. Figure 3 also shows the predicted K s from the PTFs in previous research and measured K s in this study. Figure 3 shows the range of the logarithmized K s on the left side and its probability density function on the right side.
The averaged logarithmized K s of the forest soil was 6.6, and the average logarithmized K s from the six PTFs using the same input data set varied from −2.1 to 3.6. The K S of the forest soil that was observed in this study was about 10 to 10 3 times larger than the K s predicted with the PTFs of previous research. Additionally, most of the predictions exhibited negative ln(K s ), which means that the K s was between 0 and e cm day −1 and thus very low. In this study, however, no negative values were observed. The ANOVA test also showed that the K S of the forest soil was higher than the K s predicted by previous research for forest soils, which was statistically significant. Considering the range of the soil properties in Table 4, the range of the soil size distribution is different. In particular, the bulk density of the forest soil is lower than the bulk density of the soils used in previous research. These discrepancies indicate differences of K s , and the PTFs were difficult to use to rationally explain the K S of the forest soil in South Korea.

Explanatory Variable Selection
The most appropriate treatment of the explanatory variables was selected to estimate K s through soil characteristics and topographic features. For selecting the appropriate treatment, logarithmized (ln(x)), exponential (e x ), squared (x 2 ), and untreated treatments (x) were carried out for each variable, and the linear correlation with ln(K s ) of each treatment was compared. Table 5 shows the best results of the correlation analysis of the explanatory variables. Among the soil properties, sand and organic matter had a higher linear correlation when logarithmized, and silt and clay had a higher correlation when squared (Table 5). In terms of the topographic features, the TWI was highly correlated with logarithmized treatment and CA was highly correlated when squared after being logarithmized. The p-value of all the variables was statistically significantly correlated with the logarithmized K s (p < 0.05), and the bulk density and organic matter showed a strong linear relationship.
Water 2020, 12, x FOR PEER REVIEW 8 of 19 Figure 3. Mean values and ranges of the observed logarithmized (saturated hydraulic conductivity) of forest soil in South Korea (this study), predicted logarithmized based on PTFs presented in the literature, and kernel probability density functions of these logarithmized . Six PTFs did not focus on the forest soil, and there are different ranges of the soil properties that previous research used.
The averaged logarithmized of the forest soil was 6.6, and the average logarithmized from the six PTFs using the same input data set varied from −2.1 to 3.6. The of the forest soil that was observed in this study was about 10 to 10 3 times larger than the predicted with the PTFs of previous research. Additionally, most of the predictions exhibited negative ( ), which means that the was between 0 and cm day −1 and thus very low. In this study, however, no negative values were observed. The ANOVA test also showed that the of the forest soil was higher than the predicted by previous research for forest soils, which was statistically significant. Considering the range of the soil properties in Table 4, the range of the soil size distribution is different. In particular, the bulk density of the forest soil is lower than the bulk density of the soils used in previous research. These discrepancies indicate differences of , and the PTFs were difficult to use to rationally explain the of the forest soil in South Korea.

Explanatory Variable Selection
The most appropriate treatment of the explanatory variables was selected to estimate through soil characteristics and topographic features. For selecting the appropriate treatment, logarithmized ( ( )), exponential ( ), squared (x 2 ), and untreated treatments (x) were carried out for each variable, and the linear correlation with ( ) of each treatment was compared. Table 5 shows the best results of the correlation analysis of the explanatory variables. Among the soil properties, sand and organic matter had a higher linear correlation when logarithmized, and silt and clay had a higher correlation when squared (Table 5). In terms of the topographic features, the TWI was highly correlated with logarithmized treatment and CA was highly correlated when squared after being logarithmized. The p-value of all the variables was statistically significantly correlated with the logarithmized (p < 0.05), and the bulk density and organic matter showed a strong linear relationship.  The VIF is calculated with both soil and topographic features, and values greater than 5 were found in the sand, silt, and clay fractions. This is because the sand, silt, and clay fractions are factors that are organically related to each other, so one of these factors must be removed to solve this multicollinearity problem. Sand was logarithmized, but silt and clay were squared, so sand might have a weaker linear correlation with silt and clay. Therefore, silt was removed, since silt had higher VIF values than clay. After removing silt, the VIF values of all the variables were less than 5, which indicates that there is no multicollinearity problem.

Categories with Statistical Analysis
Forest soil is affected by various forest stand and site characteristics, in addition to the soil and topographic features given in Table 6. Therefore, we analyzed the differences in K s by forest type, bedrock, landform, and soil layer ( Table 6). The forest types were divided into three types: coniferous, broadleaved, and mixed forest. Moreover, bedrock was classified into igneous, sedimentary, and metamorphic rock. Landform was classified into five categories (Table 2), and the soil layers were classified as topsoil and subsoil. According to ANOVA analysis and Tukey's post hoc analysis, the soil located in broadleaf and mixed forest had a higher saturated hydraulic conductivity than the soil located in coniferous forest (p < 0.01). Furthermore, the soil with metamorphic rock as a bedrock had a higher saturated hydraulic conductivity than the soil with igneous and sedimentary rock (p < 0.01). By landform classification, the soil located in the lower concave segment was statistically significantly higher than that in others (p < 0.01), with no statistical differences in the other four landforms except the lower concave segment. Moreover, the topsoil had a higher K s than subsoil (p < 0.01).
The overall data set was divided into several categories because there were statistically significant differences in K s by forest type, bedrock, landform, and soil layer, as shown in Table 6. As a result of the post hoc analysis of the forest type, bedrock, landform, and soil layer, all of the four forest stands and site characteristics could be classified into two categories: relatively fast K s and relatively slow K s . The total data set, therefore, was classified into 16 categories, depending on the characteristics of the forest type, bedrock, landform, and soil layer (Table 7).

Development of the Pedo-Transfer Function
The soil data were classified into 16 categories, and multiple linear regression analysis was conducted for each category. Since most PTF studies have only been conducted with soil characteristics to date, multiple linear regression analysis was conducted twice: with soil characteristics only and with soil characteristics and topographic features. Table 8 shows the results of multiple linear regression analysis with only soil characteristics as the input data. Additionally, Table 9 shows the results of multiple linear regression analysis with topographic features in addition to soil characteristics. Tables 8  and 9 show the regression coefficients of each variable according to multiple linear regression, and full equations are available in the Supplementary Materials (Table S1 and S2). Figure 4 shows the measured and predicted K s with the developed PTFs with 16 categories. The model was developed using 70% of the total data set and validated with the remaining hold-out test data set. In the calibration phase, the coefficient of determination with only soil properties was 0.37 but increased to 0.43 when topographic features were added as input data. Moreover, the coefficient of determination values, 0.37 and 0.43, are significantly different from 0. The model performance also increased in the validation phase. The RMSLE and MLE, which show the deviation of the measured and predicted values, also decreased from 0.931 and 0.724 to 0.920 and 0.709, respectively. The model performance was improved when topographic features were taken into account in both calibration and validation. Note: Clay2 * is clay2 × 10 −4 . The regression coefficients of clay2 were too small, so the constant was multiplied to make them readable. In total, 16 categories were classified according to Table 7, and these values are the regression coefficients from multiple linear regression. Note: Clay 2 * is clay 2 × 10 −4 ; elevation* is elevation × 10 −3 ; CA 2 * is CA 2 × 10 −2 . The regression coefficients of clay 2 , elevation, and CA 2 were too small, so the constant was multiplied to make them readable. In total, 16 categories were classified according to Table 7, and these values are coefficients from multiple linear regression. 0.37 but increased to 0.43 when topographic features were added as input data. Moreover, the coefficient of determination values, 0.37 and 0.43, are significantly different from 0. The model performance also increased in the validation phase. The RMSLE and MLE, which show the deviation of the measured and predicted values, also decreased from 0.931 and 0.724 to 0.920 and 0.709, respectively. The model performance was improved when topographic features were taken into account in both calibration and validation.

Model Sensitivity
Ten variables were analyzed for sensitivity analysis, and four highly sensitive variables were selected ( Figure 5). All multipliers were not applied to variables having the upper limit. For example, for the sand fraction, since the maximum value of the sand fraction is 100%, the sensitivity analysis was performed up to multiplier 4. Sensitivity analysis showed that the bulk density and sand fraction, which are soil properties, and the slope and catchment area, which are topographic features, had the great effect on . was increased when the sand fraction and catchment area increased, and

Model Sensitivity
Ten variables were analyzed for sensitivity analysis, and four highly sensitive variables were selected ( Figure 5). All multipliers were not applied to variables having the upper limit. For example, for the sand fraction, since the maximum value of the sand fraction is 100%, the sensitivity analysis was performed up to multiplier 4. Sensitivity analysis showed that the bulk density and sand fraction, which are soil properties, and the slope and catchment area, which are topographic features, had the great effect on K s . K s was increased when the sand fraction and catchment area increased, and K s was decreased when the bulk density and slope increased. However, clay, organic matter, elevation, TWI, plan curvature, and profile curvature did not significantly affect the K s .
Water 2020, 12, x FOR PEER REVIEW 13 of 19 was decreased when the bulk density and slope increased. However, clay, organic matter, elevation, TWI, plan curvature, and profile curvature did not significantly affect the .  Figure 3 shows the use of the PTFs from previous research to estimate the of the forest soil. The measured was more than 10 times larger than the predicted by previous research. Using the PTFs of previous research can therefore lead to the underestimation of the of forest soil. Because they did not focus on the forest soil and the ranges of the soil properties used for developing the PTFs were different. Unlike the soil of cultivated land, grassland, or bare land, forest soil is not compacted and large and small pores are widely distributed. Therefore, forest soil has a relatively low bulk density. As presented in Table 4, the ranges of the soil size distribution in previous research and this study are similar. On the other hand, the bulk density of the forest soil in South Korea was smaller and beyond the range of previous research. This could be one of the reasons for the higher of forest soil in South Korea. In forests, various types of plants are growing and have a great influence on the soil hydrology [25,26]. Piaszczyk et al. [27] showed that dead trees in forests can have a great impact on the soil's physical properties. Furthermore, some studies have shown that organic matter produced from the roots of trees or understory plants in forests has a significant impact on the physical structure of forest soil [22,23,28]. Due to these characteristics, forest soil has a high permeability.  Figure 3 shows the use of the PTFs from previous research to estimate the K s of the forest soil. The measured K s was more than 10 times larger than the K s predicted by previous research. Using the PTFs of previous research can therefore lead to the underestimation of the K s of forest soil. Because they did not focus on the forest soil and the ranges of the soil properties used for developing the PTFs were different. Unlike the soil of cultivated land, grassland, or bare land, forest soil is not compacted and large and small pores are widely distributed. Therefore, forest soil has a relatively low bulk density. As presented in Table 4, the ranges of the soil size distribution in previous research and this study are similar. On the other hand, the bulk density of the forest soil in South Korea was smaller and beyond the range of previous research. This could be one of the reasons for the higher K s of forest soil in South Korea.

Different Characteristics of K S in Forest Soil
In forests, various types of plants are growing and have a great influence on the soil hydrology [25,26]. Piaszczyk et al. [27] showed that dead trees in forests can have a great impact on the soil's physical properties. Furthermore, some studies have shown that organic matter produced from the roots of trees or understory plants in forests has a significant impact on the physical structure of forest soil [22,23,28]. Due to these characteristics, forest soil has a high permeability.

K S Differences by Forest Stand, Geological, and Topographical Features
Forest soil can be affected by a number of topographical factors and the type of forest stand. As shown in Table 6, it was found that differences in the K s by forest stand, bedrock, and landform were statistically significant. First, there was a difference in the K s for different forest stands, where the K s of broadleaf and mixed forests was higher than that of coniferous forest. The type of forest can have a great influence on the soil hydrology [25]. This can be explained by the root distribution difference between coniferous and broadleaf trees. Burch [24] demonstrated that coniferous trees have straight roots in a vertical direction, while broad-leafed trees spread wide in the lateral direction and produce a large number of rootlets. In other words, in the case of broadleaf and mixed forests, it is assumed that the pores of the soil are more advanced than those of coniferous forests due to the root system of trees, and thus, broadleaf and mixed forests have higher K s . Secondly, almost all of the broadleaf trees in South Korea are deciduous trees, which can produce a thicker litter layer than coniferous trees. Litter layers can produce higher organic carbon contents, which are related to lower bulk densities [27,29], so the K s increases.
Soil is made by weathering from the bedrock. Therefore, soil properties are highly affected by the bedrock. Plaster and Sherwood [30] showed that soil from metamorphic rock can include a high rate of the sand fraction. On the other hand, soil with fine particles is produced from sedimentary rock. Table 6 shows that the soil from sedimentary rock had the lowest K S , and the soil from metamorphic rock had the highest K s . The 1456 samples of forest soil also showed that the sand fracture contents were 40.6%, 35.6%, and 48.2% in igneous, sedimentary, and metamorphic rock, respectively, from which it was confirmed that soil from metamorphic rock has the highest permeability. This is consistent with the analysis of previous research [30].
The K s of forest soil also differed, depending on the landform. The lower concave segment is the part where the soil does not erode but is deposited, unlike other landforms. As the soil is deposited, the soil is not compacted, so there is a lot of space for water to move in the soil, which is the reason for the high K s . The 1456 samples showed a lower bulk density in the lower concave segment than other landforms. Martin [31] also identified that the soil property changes were those associated with the downward movement of water and soil, and there was a deeper topsoil depth and higher organic matter content in the lower concave segment, which is the lower part of the mountain. These differences can lead to high K s . Topsoil and subsoil also showed a difference in K s . Subsoil is located below topsoil, so the soil compaction derived from the gravity leads to lower K s . The subsoil soil layer also has lower organic matter contents than topsoil, which can lead to higher bulk densities and lower and lower K s [27,29]. In this study, topsoil contained 2.5% more organic matter than subsoil.
Forest soil is affected by various topographical factors, such as the forest type, bedrocks, and the landform directly related to the erosion and sedimentation. Therefore, in order to understand the characteristics of the K s , all of these factors must be considered and analyzed.

Relationship between K S and Soil and Topographic Features
The relationship between K s and soil and topographic features can be found in Tables 8 and 9. If the coefficient is positive, the variable is positively correlated with the K s . Therefore, K s increases as the sand fraction and upper catchment area increase, and K s decreases as the bulk density, elevation, and slope increase. The sand fraction and bulk density are directly related to the porosity. Therefore, a higher sand fraction and lower bulk density result in a high porosity, which can lead to a higher K s . The lower the elevation, the higher the K s , and this can be explained by the erosion and deposition process in the mountain slope. In other words, a low elevation area is also the part where the forest soil is deposited, so the soil layer is getting deeper and a higher porosity can be formed [32]. However, it can be highly related to the relationship between the elevation and other topographic features. In the lower concave segment, the K s is higher than that of other landform types (Table 6). This could be the major reason for the negative relationship between the K s and the elevation, and further research is needed to clarify the effects of the elevation on the K s . The slope and catchment area can also be explained by similar reasons. The erosion progresses as the slope increases; on the other hand, deposition progresses as the slope decreases. Additionally, as the upper catchment area increases, more water and soil are transported down the hill. In other words, the lower the elevation and the lower the slope, which is related to the bigger upper catchment area, the deeper the soil layer and the lower the bulk density, which leads to higher K s .
The organic matter, TWI, plan curvature, and profile curvature, however, did not significantly affect the K s . Nemes [33] explained that the organic matter in soil may increase the total porosity of the soil, thereby increasing the hydraulic conductivity, but in contrast, it can reduce the K s by retaining the water. In other words, one measure of organic matter contents alone cannot fully explain the K s since the effect of organic matter on the K s of forest soil is complex. The TWI, plan curvature, and profile curvature were also used as important factors for estimating the moisture content of soil, which is one of the most important soil properties [34]. However, it was found that these factors did not have a significant role in explaining the K s . Figure 4 shows that the model's performance increased when the topographic features were added, which can be confirmed through the sensitivity analysis ( Figure 5). Four highly sensitive variables are the bulk density, sand fraction, slope, and catchment area, which can explain the K s well. The bulk density and sand fraction are soil properties, but the slope and catchment area, which are topographic features, also highly affected the K s . In other words, it can be seen that K s is not only directly affected by soil properties but also greatly affected by topographic features. Therefore, topographic features should be considered when forest soil characteristic analyses are carried out.

Limitations and Suggestions for Future Research
In this study, non-linearity was considered by conducting the preprocessing of 10 soil and topographic characteristic variables to estimate K s . Furthermore, 16 categories were divided through statistical analysis and the linear relationship with K s was analyzed for each category. We confirmed that the soil properties and topographic features highly affect the K s of forest soil in this study. In South Korea, a soil investigation is being conducted on a national scale. With these database and topographic data, the K s can be estimated on a national scale by using the PTFs. Because a lot of hydraulic simulation models use the K s as an important input variable, these PTFs can be useful for running these models and for understanding the water cycle.
This model, however, still does not have a high performance. This is because the K s of forest soil is affected by more site and environmental factors than the above variables used in this study. Moreover, it is difficult to explain all of these nonlinear relationships with simple multiple linear regression. Therefore, it may be possible to construct a model with a higher performance when the other factors highly related to the K s are found and analyzed. Recently, machine learning has been used to develop a model that predicts the K s , in which case the performance of the model is increased [1,34]. However, the machine learning model is difficult to access for other researchers, so it is not easily available for other researchers. For these reasons, the machine learning model was not covered in this paper. Further research can use artificial intelligence and deep learning technology to check the nonlinear relationship between the K s and various topographical factors in forest soil, and then use several nonlinear models to develop PTFs with a high performance. In addition, it is expected that a more accurate prediction of K s on a national scale can be conducted when artificial intelligence and the deep learning model is developed using these forest soil data.

Conclusions
In South Korea, forest covers about two-thirds of the land area of the whole country. Due to these geographical characteristics, understanding the soil properties in forests is important for sustainable water management. The K s is one of the most important input variables for running a hydraulic simulation model. Therefore, developing PTFs for predicting the K s is necessary. In this study, data from 1456 sampling points were collected throughout South Korea, which are located in the temperate climate zone, and PTFs to predict the K s from soil and topographic features were developed. The K s of the broadleaf and mixed forests was higher than that of the coniferous forests, which is because of the root characteristic differences and the litter layer by forest type. The K s of the soil based on igneous and sedimentary rocks was higher than that of the soil based on metamorphic rocks. This is because of the formation of sandy soil during the weathering of metamorphic rocks. In addition, soil located in the lower concave segment had a lower bulk density and higher K s because the soil located in the lower concave segment is deposited rather than eroded, unlike other landforms. Furthermore, the K s of subsoil was lower than that of topsoil due to soil compaction by gravity and organic matter contents. Many previous researchers have only considered soil properties to predict the K s . However, since forests display spatial variability, topographic features should be considered, in addition to soil properties. The PTF model performance increased when topographic features were added. As a result, the K s of the forest soil was increased when the bulk density, clay, elevation, and slope were decreased and the sand fraction and upper catchment area were increased. The organic matter, TWI, plan curvature, and profile curvature did not highly affect the K s of the forest soil. According to the sensitivity analysis, the bulk density, sand fraction, slope, and catchment area highly affected the K s . Therefore, the topographic features were as important in predicting the K s as the soil properties of the forest soil.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/8/2217/s1, Equations S1-S6: Six PTFs were selected (Table 3) and their equation were showed. Table S1: Regression equations of categorical PTFs for the logarithmized saturated hydraulic conductivity with only soil properties. Table S2: Regression equations of categorical PTFs for the logarithmized saturated hydraulic conductivity with topographic features in addition to soil properties.