Impacts of Soil Properties, Topography, and Environmental Features on Soil Water Holding Capacities (SWHCs) and Their Interrelationships

Abstract

Soil water holding capacities (SWHCs) are among the most important factors for understanding the water cycle in forested catchments because they control available plant water that supports evapotranspiration. The direct determination of SWHCs, however, is time consuming and expensive, so many pedotransfer functions (PTFs) and digital soil mapping (DSM) models have been developed for predicting SWHCs. Thus, it is important to select the correct soil properties, topographies, and environmental features when developing a prediction model, as well as to understand the interrelationships among variables. In this study, we collected soil samples at 971 forest sites and developed PTF and DSM models for predicting three kinds of SWHCs: saturated water content (${\theta }_S$) and water content at pF1.8 and pF2.7 (${\theta }_{1.8}$ and ${\theta }_{2.7}$). Important explanatory variables for SWHC prediction were selected from two variable importance analyses. Correlation matrix and sensitivity analysis based on the developed models showed that, as the matric suction changed, the soil physical and chemical properties that influence the SWHCs changed, i.e., soil structure rather than soil particle distribution at ${\theta }_S$, coarse soil particles at ${\theta }_{1.8}$, and finer soil particle at ${\theta }_{2.7}$. In addition, organic matter had a considerable influence on all SWHCs. Among the topographic features, elevation was the most influential, and it was closely related to the geological variability of bedrock and soil properties. Aspect was highly related to vegetation, confirming that it was an important variable for DSM modeling. Information about important variables and their interrelationship can be used to strengthen PTFs and DSM models for future research.

1. Introduction

The soil water holding capacity (SWHC) is the amount of water content in soil at the particular matric suction. This is one of the most important factors for understanding and modeling the water cycle in forest catchment [1,2,3]. SWHCs are usually represented through the soil water retention curve, which is the relation curve between soil water content and applied matric suction [4]. Since this matric suction indicates gravitational force, capillary retention, and root pressure, identifying the water content of the soil corresponding to the particular matric suction is essential for understanding the water cycle, especially in relation to soil. In general, the logarithm of the absolute value of matric suction (pF) is widely used to express particular matric suction [5,6]: $\mathrm{pF}={\log }_{10}\left|\mathrm{F}\right|$, where F is the height of the water column (cm).

Water content at a specific pF value plays a major role in representing the characteristics of forest soil. Saturated water content (${\theta }_S$; water content at pF0) refers to the amount of water when the soil is saturated with water and indicates porosity, as well as the total amount of pores combined with macropores and micropores in the soil. pF1.8 indicates the starting point of capillary water holding in soil, and pF2.7 expresses the end point of the gravitation drainage [6]. pF2.7 also represents the boundary condition of easily available water for some tree species (Cryptomeria japonica and Chamaecyparis obtusa) [6]. Specifically, pF1.8 and pF2.7 are often used in forest soils, because the difference between the water content at pF1.8 (${\theta }_{1.8}$) and the water content at pF2.7 (${\theta }_{2.7}$) is directly related to the coarse capillary pore and plant growth [7]. Cianfrani et al. [5] and Wessolek et al. [8] used ${\theta }_{1.8}$ and ${\theta }_{2.7}$ to indicate the plant available soil water (PAW) and confirmed that these values were closely related to plant water uptake and evapotranspiration in the catchment scale. In addition, ${\theta }_{2.7}$ is similar to the field capacity (${\theta }_{2.5}$; water content at pF2.5 of −33 kPa), which is commonly used in soil research.

Despite its importance, estimating SWHCs is time consuming and expensive [4]. One of the most famous models for predicting SWHCs is the use of pedotransfer functions (PTFs). This is an in situ measurement-based model that predicts SWHCs based on soil property data [9]. PTF directly predicts SWHCs using soil particle distribution, bulk density, organic matter, etc. [10,11], or estimates the parameters of the soil water retention curve, indicating the relationship between the matric suction and the soil water content [12,13].

The other way to predict SWHCs is digital soil mapping (DSM). First, environmental covariates that geographically reference information corresponding to soil sampling sites are collected [14]. Prediction models are developed through environmental covariates and SWHCs, and they can be applied to large areas [15,16,17]. Topographical features, geology, and historical soil information are generally used as environmental covariates [18,19].

In order to effectively develop the PTFs and DSM models, it is important to accumulate an efficient dataset. Due to recent developments in computing technology, machine learning techniques have become popular. A supervised machine learning model is a data-driven method, and a primary way to increase model performance is to secure many training data related to the model. Thus, when various soil surveys are conducted, effective prediction model development becomes possible.

In South Korea, a nationwide forest soil investigation has been conducted by the National Institute of Forest Science (NIFoS), and SWHC prediction PTFs and DSM models were developed based on this database. Additionally, the Korea Forest Service is planning to conduct an additional nationwide soil survey to strengthen PTF and DSM models to accurately predict SWHCs in forest soil. For efficient additional data accumulation, soil surveys should be conducted considering factors that are highly influential in SWHCs. For example, if one variable (e.g., elevation or aspect, etc.) acts as an important variable in model prediction, there may be several methods that can be used. We could make that variable’s distribution of sampling sites uniform or similar to the distribution of that variable of a target area. Therefore, it is essential to select the important variables and comprehensively understand the interrelationships between soil properties, topography, environmental features, and SWHCs before conducting an additional soil survey [15]. However, there are only a few studies concerning the important variables for predicting ${\theta }_S$, ${\theta }_{1.8}$, and ${\theta }_{2.7}$ and their interrelationships in temperate forest regions.

In this context, our objectives are as follows: (1) to suggest the highly effective soil physical and chemical properties, topographies, and environmental features for the development of PTFs and DSM models and (2) to identify the interrelationships and non-linear effects on SWHC changes through selected variables.

2. Materials and Methods

2.1. Study Sites

This study was conducted in South Korea, which is the southern half of the Korean Peninsula located in East Asia. The north side of South Korea is bordered by North Korea, and the other three sides are surrounded by the sea. It has a temperate climate and four seasons. The annual average temperature is about 12.8 degrees Celsius, and the annual average precipitation is about 1343 mm. After the Korean War ended in 1953, parts of the country were destroyed, and many trees were planted as a government policy. Now, almost two-thirds of South Korea, 62.8%, is covered by forests, in a context of complex topography.

2.2. Forest Stand and Soil Properties

We used 11,544 soil samples collected at depths of 10 and 30 cm in a forested area of South Korea. These forest soil samples were collected by the National Institute of Forest Science to develop the SWHCs prediction models; soil samples were collected from 971 sample sites that were randomly selected (Figure 1). Nine hundred and seventy-one sample sites contain twelve soil texture classes, three forest types, and three kinds of bedrock (Figures S2 and S3). Six soil samples were collected in a sample site at the same point as replicates and collected at the equivalent soil depth (10 or 30 cm), though separated by a few centimeters. The soil texture distribution of the collected soil samples is shown in Figure 2.

All soil samples were collected below the organic O horizon. Since the average soil depth of A horizons in South Korean forests is about 17 to 20 cm, soil samples 10 cm deep represent mineral A horizons, while soil samples 30 cm deep represent B or C horizons.

In this study, three kinds of SWHCs are used: saturated water contents (${\theta }_S$), water content at pF1.8 (${\theta }_{1.8}$), and water content at pF2.7 (${\theta }_{2.7}$). A pressure plate method was used to measure SWHCs. The collected soil samples were analyzed by means of laboratory methods to determine six physical and chemical properties, including bulk density and organic matter. Bulk density is related to soil compaction and structure. It was measured by the ratio of the dry weight of soil and the volume of the soil. Soil organic carbon was measured by the Walkley-Black method [20], and organic matter was recalculated by multiplying by 1.724, the van Bemmelen factor. Hydraulic conductivity was measured by the falling head method, which uses a constantly changing pressure head. Soil particle fractions were measured with the hydrometer method, and three particle-size classes were classified: sand (0.05 to 2 mm), silt (0.002 to 0.05 mm), and clay (<0.002 mm). In situ measurements were also investigated to develop PTFs for analysis of the environmental impact on soil. Dominant tree height, dominant tree DBH, average DBH, and tree density were measured by forestry professionals when soil samples were collected.

Soil type was not used as an explanatory variable in this study, because of its generality and diversity issues. The soil classification system for forest soil in South Korea is totally different from other classification systems such as World Reference Base for Soil Resources (WRB) or USDA soil taxonomy. In addition, because of the relatively small territory and the similar climatic zone throughout the whole country, over 86% of the forest soils are grouped in an order, Brown Forest Soils, which is similar to Inceptisols from USDA soil taxonomy. For example, about 87.4% of forest soils used in this study are Brown Forest Soils. Because of these reasons, soil type cannot be a useful variable for predicting SWHCs in South Korea, and we did not include the soil type as an input variable.

In

Table 1. Measured forest stand and soil physical and chemical properties.

Forest Stand and Soil Physical and Chemical Properties	Abb.	Unit	At 10 cm Soil Depth (n = 953)	At 30 cm Soil Depth (n = 971)
Saturated SWC	${\theta }_S$	%	$60.4 \pm$7.9	$57.7 \pm$7.5
SWC at pF1.8	${\theta }_{1.8}$	%	$32.2 \pm$6.8	$33.0 \pm$7.2
SWC at pF2.7	${\theta }_{2.7}$	%	$26.3 \pm$5.8	$27.4 \pm$6.5
Bulk density	${\rho }_b$	g cm−3	$0.95 \pm$0.2	$1.05 \pm$0.2
Organic matter	OM	%	$4.07 \pm$1.88	$2.98 \pm$1.52
Hydraulic conductivity	${\mathrm{K}}_{\mathrm{S}}$	cm s−1	$0.015 \pm$0.013	$0.011 \pm$0.010
Sand fraction	Sand	%	$39.5 \pm$16.3	$38.4 \pm$17.2
Silt fraction	Silt	%	$37.6 \pm$15.7	$36.2 \pm$16.7
Clay fraction	Clay	%	$22.8 \pm$10.3	$25.4 \pm$11.8
Dominant tree height	DTH	m	$14.7 \pm$3.5	$14.6 \pm$3.4
Dominant tree DBH	DTD	cm	$30.6 \pm$10.2	$30.4 \pm$9.9
Average DBH	AD	cm	$24.0 \pm$7.5	$23.7 \pm$7.4
Tree density	TD	trees ha−1	$578 \pm$285	$580 \pm$283

Note: Average$\pm$standard deviation. DBH is diameter at breast height.

, soil physical and chemical properties and forest stand characteristics are shown. Some soil samples at 10 cm depth from 18 sample sites were lost. Therefore, the number of analyzed soil sampling sites is different (n = 953 for 10 cm soil depth, n = 971 for 30 cm soil depth), and the averaged values of in situ measurements at 10 and 30 cm soil depth are different (Table 1).

2.3. Environmental Covariates

We extracted an environmental covariates dataset at corresponding soil sample sites from four geographical maps: digital elevation map (DEM; 10 m resolution), geologic map (GM; 1:50,000), forest type map (FTM; 1:25,000), and forest site and soil map (FSSM; 1:25,000). Eight topographic variables from the DEM, one bedrock variable from the GM, four forest stand characteristics from the FTM, and eight forest site and soil properties from the FSSM were selected to develop the DSM models (Table S1). We chose these 4 maps because they are open source and easily available at the National Geographic Information Institute. While all the variables in the PTF dataset were continuous, the environmental covariates dataset included 11 discrete variables (Table S1).

2.4. Random Forest Model

Random forest is one of the most popular machine learning methods. Unlike other machine learning models, the random forest model is appropriate for model development, as it has discrete input variables and is based on a decision tree algorithm. There are 11 discrete variables for developing DSM models, and this is the main reason why we used a random forest model to predict SWHCs in this study. The random forest model was established using Python (v. 3.7.4) and the RandomForestRegressor module from scikit-learn (v. 0.23.2).

2.5. Variable Importance Measurements

In this study, we selected two methods for calculating variable importance: feature importance and permutation importance. Feature importance is one of the most widely used methods for variable importance measurement in the random forest model. It is calculated based on mean decrease impurity. It is related to the random forest model structure and its developing process, where the decision tree is extended by impurity [21]. In other words, the greater the decrease in impurity, the higher the importance of the variable. We calculated the feature importance using scikit-learn version 0.23.2. Permutation importance is determined by the performance differences between a basic model and a modified model. The modified version uses an input variable in which one variable is randomly permutated and the others are fixed. If a critical variable is employed to develop the modified model, the performance of the new model will decrease more than that of a less important variable. In this paper, we used the coefficient of determination to measure the performance, with 500 times random analysis on each variable.

2.6. Linear Relationship

Pearson’s correlation coefficient is used to confirm the linear relationship between explanatory variables and SWHCs. It is a method used primarily in descriptive statistics, and it has a range from −1 to 1. As the absolute value of the correlation coefficient becomes larger, the linear relationship between the two variables becomes greater. A negative value means anti-correlation and a positive value means a positive relationship between two variables. The equation for Pearson’s correlation coefficient in this study is as follows:

$$r_{XY}=\frac{{{\displaystyle \sum }}_{i=1}^N\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left(X_i-\overline{X}\right)}^2}·\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left(Y_i-\overline{Y}\right)}^2}}$$

(1)

where N is the sample size of variables X and Y.

2.7. Developed Models

A total of 12 random forest models were developed with different explanatory variables, soil sample layers, and matric suctions of response variables (

Table 2. Developed models.

Model ID	Explanatory Variables	Soil Sample Layer	Matric Suction of Response Variable
PTF-10-pF0	Forest stand and soil physical and chemical properties	10 cm depth	pF0 (saturated)
PTF-10-pF1.8			pF1.8
PTF-10-pF2.7			pF2.7
PTF-30-pF0		30 cm depth	pF0 (saturated)
PTF-30-pF1.8			pF1.8
PTF-30-pF2.7			pF2.7
DSM-10-pF0	Environmental covariates	10 cm depth	pF0 (saturated)
DSM-10-pF1.8			pF1.8
DSM-10-pF2.7			pF2.7
DSM-30-pF0		30 cm depth	pF0 (saturated)
DSM-30-pF1.8			pF1.8
DSM-30-pF2.7			pF2.7

). Explanatory variables are divided into two types: forest stand and soil physical and chemical properties for the PTF model and environmental covariates (Table S1) for the DSM model. We also classified soil sample layers at depths of 10 and 30 cm and soil water contents at saturated, pF1.8, and pF2.7 as the response variables.

2.8. Sensitivity Analysis

Sensitivity analysis can determine how much response variables are affected by explanatory variable changes. In this study, we conducted a sensitivity analysis using the developed random forest models. First, a soil sample was randomly chosen from the dataset. The explanatory variable became the input variable by changing equal intervals, while the others were kept constant. We used the modified variable as an input variable to test the sensitivity of the explanatory variable on the corresponding response variable. In this study, the analysis of each variable was repeated 150 times.

While scatter plotting and correlation coefficient analysis allow for an acknowledgement of the relationships between the explanatory variable and response variable, these analyses contain all the effects of variables that might affect response variable. On the other hand, sensitivity analysis enables the rejection of impacts of other variables on various models. Therefore, we can analyze the effects of one explanatory variable in greater detail.

3. Results

3.1. Variable Importance for Predicting SWHCs

Two methods were used for measuring variable importance. We standardized the variable importance of each model to have the total sum of 1. Feature importance and permutation importance were calculated by means of 12 models. Variable importance of forest stand, soil physical and chemical properties, and environmental covariates is shown in Tables S2–S4. To rank the importance of each variable, we averaged standardized feature importance and permutation importance (Figure 3). We also selected four critical variables for each model (

Table 3. Four primary important variables of the twelve developed models.

Model ID	First Important Variable	Second Important Variable	Third Important Variable	Fourth Important Variable
PTF-10-pF0	${\rho }_b$	OM	Sand	${\mathrm{K}}_{\mathrm{S}}$
PTF-10-pF1.8	Sand	OM	${\mathrm{K}}_{\mathrm{S}}$	DTH
PTF-10-pF2.7	Sand	OM	${\mathrm{K}}_{\mathrm{S}}$	Clay
PTF-30-pF0	${\rho }_b$	OM	Sand	Clay
PTF-30-pF1.8	Sand	${\mathrm{K}}_{\mathrm{S}}$	OM	DTH
PTF-30-pF2.7	Sand	${\mathrm{K}}_{\mathrm{S}}$	OM	Clay
DSM-10-pF0	Elevation	Aspect	STC	TSD
DSM-10-pF1.8	STC	Elevation	TPI	Bedrock
DSM-10-pF2.7	STC	Elevation	TPI	Aspect
DSM-30-pF0	Elevation	STC	Aspect	TPI
DSM-30-pF1.8	STC	Elevation	TPI	Bedrock
DSM-30-pF2.7	STC	Elevation	Bedrock	TPI

).

As illustrated in Figure 3, variable importance to predict ${\theta }_S$, ${\theta }_{1.8}$, and ${\theta }_{2.7}$, respectively, has different characteristics. In the PTFs models that used forest stand and soil physical and chemical properties, bulk density greatly impacted on ${\theta }_S$ prediction at both 10 and 30 cm. In predictions for ${\theta }_{1.8}$ and ${\theta }_{2.7}$, however, the sand fraction was the most significant. Organic matter was critical for all cases, and hydraulic conductivity (K_S) was also important. In the PTF-10-pF2.7 and PTF-30-pF2.7 models, the influence of clay increased, and the ranking of dominant tree height (DTH) on PTF-10-pF1.8 and PTF-30-pF2.7 went up, while the variable importance value remained low (below 6%).

In the DSM models that used environmental covariates as input variables, elevation was the primary important variable in ${\theta }_S$ prediction at soil depths of 10 and 30 cm. To predict ${\theta }_{1.8}$ and ${\theta }_{2.7}$, the soil texture class was the most significant. The next important variables were TPI and bedrock. In the DSM-10-pF0 model, aspect showed a high variable importance.

From the variable importance analysis, the importance difference by soil depth layer (10 and 30 cm) was not significant; however, the difference shown by different matric suctions was larger. The order of important variables between ${\theta }_S$ and ${\theta }_{1.8}$ showed a large difference, and the order of important difference between ${\theta }_{1.8}$ and ${\theta }_{2.7}$ was small, showing similar important variables (Table 3). The forest stand dataset influenced less than 11.6% in the PTF models (average 7.6%); therefore, soil physical and chemical properties were dominant in these models. In the DSM models, the digital elevation model (DEM) and forest site and soil maps (FSSM) showed more than 80% importance.

3.2. Correlation between Highly Effective Variables and SWHCs

From the variable importance analysis, we selected eight highly effective explanatory variables: four forest stand and soil physical and chemical properties (${\rho }_b$, K_S, sand, and OM) and four environmental covariates (elevation, STC, TPI, and bedrock). To identify the linear interrelationship between highly effective input variables and SWHCs, we plotted correlation matrix plots (Figure 4).

First, three kinds of SWHC were closely interrelated. The correlation coefficients of each relationship among the SWHCs were all positive, especially ${\theta }_{1.8}$ and ${\theta }_{2.7}$. In soil physical and chemical properties, bulk density (${\rho }_b$) and sand fraction (sand) had a negative influence on SWHCs. Bulk density and organic matter showed a negative correlation with a high coefficient, and organic matter was positively related with SWHCs. All other variables except K_S presented the same tendency to the three kinds of SWHCs. K_S showed a positive correlation with ${\theta }_S$, but a negative correlation with ${\theta }_{1.8}$ and ${\theta }_{2.7}$.

Regarding environmental covariate, elevation and TPI presented a positive relationship with SWHCs, and soil texture class and bedrock, which have a high correlation coefficient value to sand, showed a negative relationship with SWHCs. Elevation showed a high correlation with three soil properties. It had a negative relationship with bulk density and sand fraction and a positive relationship with organic matter. The differences in interrelationships between selected variables and SWHCs at soil depths of 10 cm and those at soil depths of 30 cm were not significant.

3.3. Sensitivity Analysis for Identifying Non-Linear Relationship

To identify the non-linear effect of input variables on SWHCs, a sensitivity analysis was conducted. We selected four variables that were primary and secondary important variables in PTF and DSM, respectively, and analyzed sensitivity on ${\theta }_S$ and ${\theta }_{2.7}$ (Figure 4). Since ${\theta }_{1.8}$ and ${\theta }_{2.7}$ are closely related (see in Figure 4) and have a similar trend, we excluded ${\theta }_{1.8}$ from Figure 5.

In Figure 5a–h, ${\theta }_S$ was highly sensitive to bulk density and elevation changes. As shown in Figure 4, bulk density has a negative relationship with ${\theta }_S$, and a positive relationship with elevation. Notably, ${\theta }_S$ steeply increased at the elevation between 500 and 1000 m. On the other hand, sand fraction and soil texture classes showed a modest contribution when predicting ${\theta }_S$ and presented a low sensitivity.

However, ${\theta }_{2.7}$ showed the opposite tendency. Bulk density and elevation, which had high sensitivity in ${\theta }_S$, showed low sensitivity in ${\theta }_{2.7}$, and sand fraction and STC showed high sensitivity.

Figures S3–S6 depict sensitivity analysis results of the four selected variables from twelve developed models. From the sensitivity analysis, we confirmed the non-linear effects on SWHCs, not merely the linear relationship between two variables. For example, organic matter was significant only when it was below 5%. Hydraulic conductivity also showed a significant relationship only when K_S < 0.02 cm s⁻¹. Moreover, the impact of aspect on saturated soil water content was limited at a depth of 10 cm. The southern aspect (90$°$ to 270$°$) had a lower ${\theta }_S$ than the northern aspect (0$°$ to 90$°$ and 270$°$ to 360$°$), the tendency of which was hardly detected by a linear relationship such as Pearson’s correlation.

4. Discussion

4.1. Influential Soil Physical and Chemical Properties on SWHCs Prediction

The primary important soil property to predict ${\theta }_S$ was bulk density, while sand fraction was the most valuable in predicting ${\theta }_{1.8}$ and ${\theta }_{2.7}$ (Figure 3). Sensitivity analysis also demonstrated that the impact of bulk density and sand fraction on soil was different, along with different matric suction (Figures S3 and S4). Since bulk density is the mass of soil per unit volume, it represents soil structure and soil compaction [22]. The result of sensitivity analysis showed that an increase in the bulk density led to a decrease in ${\theta }_S$. Thus, compacted soil has fewer macropores and a lower water holding capacity.

On the other hand, sand fraction was the first important variable for predicting ${\theta }_{1.8}$ and ${\theta }_{2.7}$, not silt or clay, which both showed a negative relationship. Since ${\theta }_{2.7}$ is the amount of soil moisture after being pressured from pF1.8 and pF2.7 (−6.2 kPa and −49.1 kPa, respectively), soil water in macropore is eliminated at ${\theta }_{2.7}$ [23]. Therefore, the trivial influence of bulk density on ${\theta }_{1.8}$ and ${\theta }_{2.7}$ prediction (below 4%) is explainable by soil texture, which is related to capillary force and micropore rather than macropore and soil structure.

The PTF-10-pF2.7 and PTF-30-pF2.7 models denoted that importance of the sand fraction became greater than that in predicting ${\theta }_{1.8}$, and it was also shown in sensitivity analysis results. Clay was significant in ${\theta }_{2.7}$ prediction. Water holding capacity is closely related to capillary force [23]. As the particle of soil becomes finer, its surface area becomes broader; therefore, the capillary is enhanced. Soil needs a wider surface area and finer soil particles such as silt and clay to have more moisture at high matric suctions. In this regard, the enhanced role of the sand fraction in predicting ${\theta }_{2.7}$ is indirect. The sand fraction is equal to the summation of the silt fraction and the clay fraction in hydrological modeling, since the sum of sand, silt, and clay is 100%. Small particles such as silt and clay become powerful in ${\theta }_{1.8}$ and ${\theta }_{2.7}$ predictions. Moreover, a higher rank of variable importance with a weight of 5~7% emphasizes the influence of clay in ${\theta }_{2.7}$, whereas the variable importance of the clay fraction for predicting ${\theta }_S$ and ${\theta }_{1.8}$ is minute (Figures S3 and S4).

In this vein, in order to estimate water holding capacity at a higher matric suction, such as ${\theta }_{4.2}$ (water content at pF4.2; permanent wilting point), the importance of the clay fraction increases. Qiao et al. [24] employed a stepwise method to select proper variables for predicting ${\theta }_{2.5}$ (field capacity) and ${\theta }_{4.2}$ (permanent wilting point). While the sand fraction was used as the main factor in the ${\theta }_{2.5}$ prediction, they selected clay as a main factor in the ${\theta }_{4.2}$ prediction. Adhikary et al. [25] also used the clay fraction for ${\theta }_{4.2}$ prediction, whereas the sand and silt fractions were used for ${\theta }_{2.5}$. As the capillary force is closely related to the soil texture (bulk density < sand < silt < clay), corresponding soil particle distribution mainly affects the SWHCs (${\theta }_S$ < ${\theta }_{1.8}$ < ${\theta }_{2.7}$ < ${\theta }_{4.2}$). In other words, soil structure is the main factor in ${\theta }_S$ rather than soil texture, sand and silt are the major factors in ${\theta }_{1.8}$, and silt and clay are the most influential variables in ${\theta }_{2.7}$.

The relationship between the K_S and SWHCs can be explained by the relation mentioned above. In Figure 4, K_S has a positive correlation with ${\theta }_S$, but a negative correlation with ${\theta }_{1.8}$ and ${\theta }_{2.7}$. K_S and bulk density are negatively correlated, and compacted soil has smaller K_S, since it has fewer macropores. On the other hand, the bulk density effect is negligible and soil texture is dominant in ${\theta }_{1.8}$ and ${\theta }_{2.7}$. Therefore, a higher value of K_S can be explainable through a high sand fraction content or a small clay fraction content, which are both closely related to the smaller value of ${\theta }_{1.8}$ and ${\theta }_{2.7}$. The scatter plot in Figure 4 shows a weak correlation between K_S and sand fraction; this is because K_S is mainly affected by disturbances in the soil structure, indirectly caused by factors that affect bulk density and organic matter [22].

Organic matter was another powerful factor in the prediction of SWHCs. In Figures S3 and S4, SWHCs rise as organic matter contents become greater. The influence of organic matter on ${\theta }_S$ is linked to bulk density. Organic matter enriches porosity and lessens soil compaction [26,27,28,29]. In particular, the negative relationship between organic matter and bulk density has been suggested as a linear equation in previous studies [30]. We also confirmed the linear relationship between bulk density and organic matter (BD = 1.26–0.073 OM; Figure 4). In other words, an increase in organic matter results in a reduction in bulk density and expansion of macropore; therefore, ${\theta }_S$ rises.

Organic matter also affects water holding capacities. Many studies have confirmed that organic matter can increase ${\theta }_{2.5}$ [3,31,32,33], which is closely related to ${\theta }_{1.8}$ and ${\theta }_{2.7}$ [5]. Leu et al. [32] confirmed that organic matter even increases in ${\theta }_{4.2}$. SWHCs are also highly related to the plant available water (PAW; generally represented as ${\theta }_{2.5}$-${\theta }_{4.2}$ or ${\theta }_{1.8}$-${\theta }_{2.7}$). Lal [34] demonstrated that the influence of organic matter on capillary force was more noticeable at field capacity (${\theta }_{2.5}$) than at permanent wilting point (${\theta }_{4.2}$), raising PAW. In this study, we also found similar results in our sensitivity analysis. In Figures S3 and S4, sensitivity on organic matter is higher at ${\theta }_{1.8}$ than ${\theta }_{2.7}$. Thus, we identified that organic matter was effective in lower matric suction. We also confirmed that, when organic matter content is more than 5%, sensitivity dropped drastically. It is thought that organic matter with a large particle size is too large to affect capillary force. Consequently, soil organic matter shrinks bulk density and raises ${\theta }_S$ and capillary force, as well as ${\theta }_{1.8}$ and ${\theta }_{2.7}$. In this regard, we found that organic matter influenced SWHCs and PAW significantly. Its sensitivity was active below 5%, and it decreased sharply above 5%.

4.2. Interrelationship between Topography, Soil Properties, and Vegetation

Topography was also an influential factor in SWHCs. In previous studies, elevation, which was included in DEM, was the most dominant factor in predicting SWHCs [15,35,36]. In Figure 4, we confirm the interrelationships between elevation and soil physical and chemical properties. Elevation was closely related to bulk density, sand fraction, and organic matter of the soil. The sand fractions decreased and organic matter contents increased as the elevation became high. Most of Baekdudaegan, which is a mountain spine stretching from north to south in South Korea [37], is a nationally protected area. In this region, development is strictly limited, and natural vegetation is widely spread. It is considered that soil in high elevation has low bulk density and high organic matter contents because of these reasons. In addition, the effects of elevation on SWHCs can be found by the spatial distribution of bedrock. Since soil particles were made from bedrock, soil texture was highly related to its bedrock. Plaster et al. [38] showed the differences in soil textures from different bedrocks and sedimentary and metamorphic rocks. Metamorphic rock tends to create sand fraction, while sedimentary bedrock tends to create finer fractions. Soil samples in this study showed that the average elevation of sample sites with metamorphic rock was 286.7 m, and the average elevation of sample sites with sedimentary rock was 538.9 m. The average elevation of two types of bedrock is different and statistically significant (p < 0.01). This shows that finer fractions are more distributed in higher elevations due to the large distribution of sedimentary bedrock at high elevations.

In Figures S5 and S6, aspect was the most important topographical feature in the DSM-10-pF0 and DSM-30-pF0 models. Sensitivity analysis showed that the ${\theta }_S$ of north-facing slopes was greater than the ${\theta }_S$ of south-facing slopes. Solar energy in south-facing slopes, which is larger than that of north-facing slopes, raises soil and air temperatures of the mountainside. This could create water stress in plants and hinder the growth of vegetation, since high temperature increases evapotranspiration [39]. However, in South Korea, relatively high summer rainfall and annual rainfall, which is 1343 mm, generally limit water stress, and higher solar energy in south-facing slopes helps vegetation grow well. Vegetation and soil chemical weathering are highly related, and it can make soils in south-facing slopes finer [40,41,42]. Freeze–thaw processes could be another reason. In the winter season, the soil water freeze–thaw cycle is more frequent on south-facing slopes because of higher solar energy during the day; this cycle also promotes physical soil weathering [43,44], which could contribute to a greater relative abundance of finer particles and low ${\theta }_S$ on south-facing slopes.

4.3. Limitations and Recommendations for Future Research

In this study, we collected data related to soil physical and chemical properties, environmental covariates, and forest stand characteristics to investigate the contribution factors of SWHCs. Forest stand characteristics data showed marginal effects on SWHC predictions compared to other variables, such as soil properties and topographical features [35]. It seems that the role of vegetation on soil moisture capacity is subtle; however, this is because we only used forest stand characteristics for biological factors. Understory vegetation plays a key role in affecting the moisture content of topsoil, organic matter content, and regional evapotranspiration [45]. In future research, the data of understory vegetation and NDVI, which represents overall vegetation, should be included to clarify the relationship among the soil physical and chemical properties, vegetation, and environmental covariates.

We used bulk density and soil texture for soil physical properties, and these variables showed great importance in predicting SWHCs. However, coarser soil particles over 2 mm were not considered for developing the models in this study. In future research, it may be more appropriate to include coarser soil particles to develop models and to confirm their interrelationships, since different size particles showed markedly different properties, as shown in Section 4.1.

5. Conclusions

The developed pedotransfer function (PTF) models based on forest stand and soil physical and chemical properties showed that bulk density had the greatest influence on predicting saturated water content (${\theta }_S$), while sand content had the greatest influence on predicting water content at pF 1.8 and pF2.7 (${\theta }_{1.8}$ and ${\theta }_{2.7}$). The digital soil mapping (DSM) models developed using environmental covariates as an input dataset showed that elevation was the most influential factor in predicting ${\theta }_S$, and soil texture class was the most influential factor in predicting ${\theta }_{1.8}$ and ${\theta }_{2.7}$.

Variable importance and sensitivity analysis showed that, as the matric suction changed, the soil physical and chemical properties that mainly influence the soil water holding capacities (SWHCs) changed to the following values: soil structure in ${\theta }_S$; sand and silt fraction in ${\theta }_{1.8}$; and much finer particles in ${\theta }_{2.7}$. It was confirmed that the organic material increased ${\theta }_S$ by reducing the density and also increased ${\theta }_{1.8}$ and ${\theta }_{2.7}$ by increasing capillary force. The sensitivity of SWHCs to organic matter was significant when it was less than 5%. Elevation was closely related to the geological variability of bedrock and soil physical and chemical properties, and aspect was highly related to vegetation, confirming that it was an important variable for developing the DSM model.

This study contributes to the data collection process for the development of more accurate PTF and DSM models by presenting important variables necessary for the estimation of three kinds of SWHC. Unfortunately, we could not find a significant relationship between forest stand characteristics and SWHCs. However, we were able to confirm the close association with topographic features, soil physical and chemical properties, and vegetation. More accurate SWHCs prediction models will be developed when the data related to understory, NDVI data, and coarser soil particle are added.