Predicting and mapping soil available water capacity in Korea

Geography and soil data

This study considered the mainland of South Korea, excluding some smaller islands, such as Ulleung, and Dokdo. The area is approximately 100,000 km², lies between 33.931°N and 38.705°N and 125.436°E and 129.705°E. Korea is located in the humid temperate climatic zone, affected by the influence of both continental and oceanic air masses. National mean annual rainfall is 1,300 mm, ranging from 1000 mm (Daegu area) to 1,800 mm (Jeju area). Approximately half of the annual rainfall occurs during the summer months, June to August, with occasional typhoons. During the summer, although the ambient temperature is high (mean temperature; 20–25°C), and crop canopies are thick, the precipitation exceeds the potential evapotranspiration because of concentrated heavy rainfall. As a result, the base saturation ratios in the majority of soils are rather low. Temperatures in the spring and autumn are mild (mean temperatures of 10–15°C) and winter months are rather cold (monthly mean temperature of −5–2°C), particularly in the central and northern regions. In the winter, spring and autumn, the precipitation is less than in the summer. Despite these seasonal differences, the amount of precipitation and the potential evapotranspiration remain similar year round because of the lower ambient temperatures in the drier seasons.

Korea is a mountainous country. More than two-thirds of the country is occupied by mountains with steep slopes. Plains are subdivided into inland plains, coastal plains and the plains in the narrow valleys. The plains had been under intensive use for agricultural production. The high relief of the land coupled with the heavy downfalls of rain in the summer affects the characteristics of Korean soil very profoundly. Soil erosion has been intense throughout the country for a long time, particularly where the population density is high.

The parent materials of Korean soils are part of the ten recognised geologic systems from different geological time series. Dominant rock types include granitic gneiss (32.4%), granite (22.3%), schist (10.3%), limestone of the Chosun series (Cambrian-Ordovician; 10.1%). The former three geological systems are present in about 60% of the land area and are known as acidic rocks. The fact that the rainfall exceeds the potential evapotranspiration and coupled with the abundance of acidic rocks, results in the wide occurrence of acidic soils in the country. The limestones of the Chosun series are alkaline, thus the soils derived from these rocks tend to be neutral or slightly alkaline. These soils are only found near Kangwon Province. Some areas contain sandstone bedrock which results in coarse-textured soils. However, even among the soils derived from the same parent rock, the textures can vary depending upon the location in the soil catena. Soils developed in high elevation areas tend to be coarse due to severe loss of fine particles by erosion, while the soils developed at the locations where soil erosion is not severe tend to be fine-textured. Given, the high variation of soil as function of climate, topography, parent materials and landuse, it is a challenge to produce estimates of soil properties at a fine resolution.

The soil survey was initiated by the Korean Rural Development Agency (RDA), UN, and FAO with a reconnaissance survey, making use of aerial photographs purchased from USA funded by Korea Soil Survey Organization between 1964 and 1967. As a result, soil maps of Korea at scales of 1:250,000 and 1:50,000 were published. Thereafter, RDA adopted the US Soil Taxonomy system and carried out the detailed soil survey between 1968 and 1990 alone. Now, detailed soil maps (1:25,000) are available for entire country in both hard copies and digital format. Furthermore, RDA prepared highly detailed soil maps (1:5,000) between 1995 and 1999 for the entire country. These were digitized and made available for the public through the RDA web site (http://soil.rda.go.kr).

The soil database used in this study was compiled based on the soil profiles in “Taxonomical Classification of Korean Soils” (NIAST, 2000), which were mostly collected in the 1970s. It includes soil profile descriptions along with physical and chemical analysis for the modal profiles of the 380 soil series defined in South Korea. Soil chemical and physical properties of each horizon (n = 1, 559) were recorded, including particle size distribution, moisture retention, organic matter content (OM), cation exchange capacity (CEC), and a limited number of bulk density (BD) data. Soil water retention was recorded for water content at −10, −33 and −1500 kPa on the mass basis (w in g 100 g⁻¹) using soils that were ground and sieved to <2 mm. The statistics of the data is given in the Table S1. The data show a wide range of distribution for bulk density and water retention. The low bulk density and high water retention soils are due to the Andisols, while the high bulk density and low water retention soils are due to Entisols which have a high sand content (>70%).

A couple of problems arise in using this database to map the AWC:

• Not all soil samples have measured bulk density and w at −10 and −1500 kPa. Thus they need to be predicted using pedotransfer functions

• Water content (w) at −10 kPa in the database was measured using disturbed (ground) samples. Thus a correction factor is needed to take into account the effect of bulk density.

Prediction of bulk density

Only 108 samples contain measured BD, so first we derived a linear model predicting bulk density as a function of organic matter (OM) and sand content. While this model has a reasonable predictive capacity (R² = 0.59), when we applied this model to the whole datasets, we obtained unreasonable estimates of values between −2.00 to 2.50 g cm⁻³. Therefore we used a more conceptual function as proposed by Tranter et al. (2007) which first estimated the mineral bulk density, BD_min as a function of sand content and soil depth. This estimate is then combined with the BD model of Adams (1973) which incorporated the effect of organic matter (OM). Andisols which contain a high amount of allophanic materials can have a different response to the model. Therefore we derived a separate model for BD for Andisols using a global dataset (Tempel, Batjes & van Engelen, 1996).

Pedotransfer functions predicting field capacity and wilting point

Available water capacity is defined here as the amount of water held by the soil between field capacity and wilting point. Water content at field capacity is usually measured in laboratory at a potential of −10 or −33 kPa. Here we use water content at −10 kPa to represent field capacity, and at −1500 kPa for wilting point. Although field capacity is more of a dynamic property which depends on the soil’s hydraulic conductivity (Romano, Palladino & Chirico, 2011), because of the unavailability of field measurement, we only used the laboratory data.

In the absence of laboratory measurement, water content at field capacity can be predicted using pedotransfer functions from the soil’s particle size distribution, and bulk density or soil structural information (McBratney et al., 2002; Saxton & Rawls, 2006). Here we developed a linear regression model predicting gravimetric water content at −10 kPa (w₁₀) and −1500 kPa (w₁₅₀₀) measured on ground soil samples using basic soil information: clay content, sand content, organic matter content (OM) and cation exchange capacity (CEC): ${\displaystyle w_{10,1500}=f\left(\mathrm{sand,clay,OM,CEC}\right).}$

Adjustment of field capacity according to bulk density

Water content at field capacity is affected by macroporosity and soil structure (Sharma & Uehara, 1968), and therefore measurement is recommended using natural soil clods. Meanwhile water content at wilting point or −1500 kPa is not much affected by structure, as most water is held with adsorptive forces, thus it can be measured using disturbed soil samples (Aina & Periaswamy, 1985). However in the Korean soil database, water content at −10 kPa was measured on ground samples. This is because the samples collected from soil survey were mainly for mapping and classification purposes, and samples for bulk density and soil clods were not collected. This problem is also common in other countries (Bell & Van Keulen, 1996). Bell & Van Keulen (1996) warned against the use of field capacity data derived from disturbed samples as its measurement overestimates in-situ field capacity for most soils except for coarser textured soils.

Thus water retention data in the Korean soil database which were measured on ground samples need to be adjusted to represent the likely water content at a given bulk density. To build such a model, the only publicly available data that contained such measurements is the soil characterization and profile data from the US National Soil Characterization database (Soil Survey Staff, 1997). 301 samples in the database contained water content at −10 kPa measured using both natural clods and disturbed samples. From this subset, 274 samples were selected for building the model and the rest (27 samples) were used as validation data. The data are from 141 profiles, and the samples come from A, B and C horizons from various depths in the top 2 m of the profile.

The statistics of the soil properties from the US database is given in Table S2. The database contained measurements of w_{10 clod}, the percent mass of water retained at suction of 10 kPa, which was measured on natural fabric (clods), and reported on a <2 mm base. w_{10 ground} is the gravimetric water content of air dry <2 mm (ground) samples, after equilibration at 10 kPa suction. BD is the bulk density (g cm⁻³) of the <2 mm fraction, with volume being measured after equilibration at −10 kPa. We used a linear regression to obtain the estimates of w_{10 clod} from w_{10 ground} plus other basic soil properties.

Prediction of profile available water capacity

After adjustment of water retention at −10 kPa for bulk density, available water capacity for each layer is calculated as the difference in volumetric water content (θ) between field capacity and wilting point, adjusted by gravel content: (1)${\displaystyle \text{AWC (mm)}=({\theta }_{10}-{\theta }_{1500})/100\times (1-R_v)\times \text{Thickness of layer (mm)}}$ where volumetric water content (θ in percent volume) is calculated from gravimetric water content (percent mass) multiplied by bulk density (in g cm⁻³): (2)${\displaystyle \theta =w\times \mathrm{BD}.}$R_v is the volume fraction of gravel (m³ m⁻³), calculated from Saxton & Rawls (2006): (3)${\displaystyle R_v=\frac{\alpha R_w}{1-R_w\left(1-\alpha \right)}}$ where R_v is the mass fraction of gravel (kg kg⁻¹), and a is the ratio of soil to gravel density (BD/2.65). The PAWC (profile available water capacity in mm) is calculated as the sum for all the layers up to 100 mm: (4)${\displaystyle \mathrm{PAWC}=\sum _{i=1}^n{\mathrm{AWC}}_i.}$ The procedure for estimation of PAWC is summarised in Fig. 1.

Mapping

The map of the available water capacity to 1 m depth (PAWC) was made by calculating the modal PAWC value for each of the South Korean soil series. The modal values for each soil series were then allocated to each of the map units of the 1:25,000 map. An equal-area spline was utilised to represent the continuous soil depth function of AWC. The equal-area spline function not only fits the water content data with depth, but also disaggregates data obtained from horizon bulk samples into a continuous depth distribution. The spline function consists of a set of local quadratic functions tied together with ‘knots’ that describe a smooth curve through a set of points. These functions are advantageous for harmonizing soil properties that were collected at different depths. See Malone et al. (2009) for the mathematical detail of the functions.

Map validation

An independent dataset was used to validate the map predictions. The map validation dataset comprises measurement from 170 soil profiles representing 165 soil series. These profiles were collected recently in 2009 and 2011. As the samples were collected for soil survey purposes, no water retention measurement was made, however the data contains measurement of particle size analysis (850 samples) and bulk density (570 samples). We can therefore validate the maps of bulk density, clay and sand content. These variables are the important predictors of available water capacity. The validation was performed on 5 standard depths where the soil samples were measured: 0–5, 5–15, 15–30, 30–60, and 60–100 cm. The prediction from the map was resampled to these standard depths using the equal area spline function.

To validate the prediction of water content at field capacity, another published dataset was used (Han et al., 2008a; Han et al., 2008b). This dataset, here called water retention validation dataset, is a collection of 18 mostly topsoil (0–20 cm) samples, representing 7 soil series. The samples were measured for soil hydraulic properties using an undisturbed soil core with the onestep outflow technique, and water content at −10 kPa was calculated (see Table S3 for the description of the dataset).

PTF for bulk density

Following Tranter et al. (2007), we first derived a pedotransfer function (PTF) predicting mineral bulk density (BD_min) from sand content and depth based on the limited data in the database (n = 108): (5)${\displaystyle {\mathrm{BD}}_{\min }=1.017+0.0032\ast \mathrm{Sand}+0.054\ast \mathrm{log(depth)}}$ where depth is the mean depth of the sample (cm) and Sand is the percentage by mass of sand (g 100 g⁻¹). The coefficient of determination (R²) is quite low (0.20) as this model only accounts for the influence of the mineral component and overburden. The influence of organic matter on bulk density is then calculated based on the model of Adams (1973): (6)${\displaystyle \mathrm{BD}=\frac{100}{\frac{\mathrm{OM}}{{\mathrm{BD}}_{\mathrm{OM}}}+\frac{100-\mathrm{OM}}{{\mathrm{BD}}_{\min }}}}$ where OM is the mass percentage of organic matter (g 100 g⁻¹), and BD_OM = average organic matter bulk density = 0.224 g cm⁻³, and BD_min is the mineral bulk density (g cm⁻³).

Because of the different mineralogical composition of Andisols, a separate model that predicted BD from OM content was derived from the Tempel, Batjes & van Engelen (1996) dataset: (7)${\displaystyle \mathrm{BD}=1.02-0.156log\left(\mathrm{OM}\right)}$ (R² = 0.45, RMSE = 0.26, n = 642).

The relationship between OM and bulk density for the soil in the database is given in Fig. 2, which shows bulk density is decreasing with increasing OM content. The relationship for Andisols (Eq. (7)) shows lower bulk density values compared to the other mineral soils. This is to be expected as Andisols are dominated by allophanic materials which have lower density. Applying Eqs. (5), (6) and (7) to the Korean soil database provides an R² = 0.49 (n = 108).

We validated the bulk density PTF using the independent map validation dataset. Table 1 shows the validation statistics in terms of mean error (ME), root mean squared error (RMSE), and coefficient of determination (R²). The results showed that the overall prediction error (RMSE = 0.208 g cm⁻³) is higher compared to values obtained by others. For example Martin et al. (2009) obtained an RMSE = 0.123 g cm⁻³ for soils from France and Tranter et al. (2007) obtained an RMSE = 0.195 g cm⁻³ for soils in Australia. However, our prediction is not biased (ME values close to 0). The prediction is considered reasonable in this context as the pedotransfer functions were calibrated from a limited observation (n = 108). The prediction for Andisols (Eq. (7)) appeared to be useful with no bias and low RMSE.

PTF for field capacity and wilting point

From the Korean soil database, we derived a linear model predicting water content at −10 kPa (w₁₀) to represent field capacity: (8)${\displaystyle w_{10}=33.18-0.188\mathrm{Sand}+0.918\mathrm{CEC}+3.578\mathrm{log(OM)}}$ (R² = 0.66, RMSE = 8.02 g 100 g⁻¹, n = 1203).

Similarly, we derived another model for water content at −1500 kPa (w₁₅₀₀) to represent wilting point: (9)${\displaystyle w_{1500}=3.13+0.186\mathrm{Clay}+0.541\mathrm{CEC}+1.708\mathrm{log(OM)}}$ (R² = 0.61, RMSE = 4.67 g 100 g⁻¹, n = 1435).

Note that these PTFs are based on measurements from ground samples with units in percentage mass (g 100 g⁻¹). The goodness of fit for these PTFs is shown in Supplemental Information. These PTFs are used to fill the data gap of unmeasured w₁₀ which are only 25% of the dataset, and unmeasured w₁₅₀₀ which constituted 10% of the data.

Adjustment of field capacity

Using the dataset from US National Soil Characterization database (Soil Survey Staff, 1997), a relationship between w_{10 clod} and w_{10 ground} is derived (10)${\displaystyle {\mathrm{w}}_{\text{10 clod}}=40.71+0.67w_{\text{10 ground}}-21.36\mathrm{BD}}$ (R² = 0.81, RMSE = 4.45 g 100 g⁻¹, n = 274).

This model describes the reduction of soil water content with increasing bulk density (Fig. 3). We validated this model on 27 soil samples in the US database that were not used for deriving the models. For water content at −10 kPa, Eq. (10) gives an R² = 0.73. It can be seen that our models predict very well the likely water retention at −10 kPa of clod samples given data obtained from disturbed soil samples. This equation allows the estimation of the representative water content at −10 kPa for a given bulk density.

Available water capacity

Based on Eqs. (5) to (10), volumetric water content (θ) at −10 kPa, −1500 kPa, and their difference (AWC) were calculated for each soil series. Table 2 shows the predicted values of these soil properties grouped by soil order. Histosols show the highest water retention at −10 and −1500 kPa. Andisols, which are dominated by non-crystalline allophanic minerals, show the highest water storage capacity of 29%, which is also known in the literature (Shoji et al., 1996). They also have the highest volumetric water content at field capacity (mean 56%), and the volumetric water content at wilting point (mean 27%). This pattern of high water retention is also observed by Fontes, Gonçalves & Pereira (2004) and Pochet et al. (2007). As indicated by Shoji et al. (1996), one of the special characteristics of Andisols is their high water-holding capacity at −1500 kPa. Entisols and Inceptisols which have loamy texture also show high water retention. The lowest water retention is Entisols which are dominated by sandy materials.

The predicted values were compared with measured AWC values from the US soil data in the Tempel, Batjes & van Engelen (1996) database. Table 3 shows the AWC of the US soil data grouped by mineralogy and texture. For mineral soils, the highest water capacity is found in allophanic soils with a mean of 27%, similar to our estimates (Table 2). Overall medium textured soils have higher AWC in both soil with low activity and high activity clay. These values are in accordance with our estimates for Korean soils.

Finally, PAWC (amount of available water that can be stored in the soil to 1 m) was calculated using Eq. (1) for the modal soil profiles of the 380 South Korean soil series. Table 4 shows the statistics of the estimates grouped by soil order. Note that the Mollisols and Histosols were excluded in this analysis as they are only represented by 2 soil series. The soils with the highest PAWC are the Ultisols and Alfisols (mean of 206 and 205 mm, respectively). Interestingly, Andisols appear to have a low PAWC, contradictory to the horizon measurement. The low values are due to the higher gravel content of Andisols which affected the calculation of PAWC.

PAWC map

The value of PAWC for each of the soil series was then mapped onto the 1:25,000 soil map and shown in Fig. 4. The western part of the country generally showed higher available water capacity than the eastern part of the country which is mountainous and has shallower soils. The mean value of PAWC (to 1 m depth) for Korean soil is 189 mm. Total available water capacity is summarized by agricultural land use using a map provided by the Korean Ministry of Land, Transport, and Maritime Affairs. PAWC is the highest in paddy fields, with a mean of 231 mm, which are mostly located in the fluvio-marine plain (See Supplemental Information).

We also are able to display the volumetric water content at −10 and −1500 kPa down the profile as a continuous function using the equal-area spline. Examples for several soil series are shown in Supplemental Information.

Map validation

In this paper, each mapping unit (from 1:25,000) is represented by a soil series, and each soil series is represented by a modal profile. Obviously there is variability within a mapping unit, i.e. the soil is not homogenous. Here each soil series is only represented by a modal profile and we do not have an uncertainty measurement of the variability. In order to gauge the accuracy of the map, we compared the map prediction of bulk density, sand and clay content (which were used as inputs for prediction of AWC) with measured values from an independent dataset. Tables 5–7 shows the accuracy measure in terms of mean error (ME), root mean squared error (RMSE), and coefficient of determination (R²).

The map of bulk density (Table 5) is realistic with an overall R² value of 0.47, and has the same magnitude of error compared to the pedotransfer prediction (Table 1). There is an indication that the prediction at the surface layer (0–5) cm is less accurate, as the surface layer is mostly affected by cultivation. Table 6 shows that the map predicts sand content very well with R² values between 0.5 and 0.6 for all depths, meanwhile prediction of clay content (Table 7) is less accurate with R² values between 0.2 and 0.5.

Leenhardt et al. (1994) evaluated the accuracy of soil maps in an area of 1328 ha, the south of France. They showed that the R² values, which were calculated as the proportions of variance explained by the 1: 25,000 soil map, for w₁₀, w₁₅₀₀ bulk density, clay, and sand content are 0.77, 0.71, 0.73, 0.62, and 0.66, respectively. In comparison, our validation shows lower R² values for bulk density, clay, and sand content (R² values of 0.47, 0.56, and 0.40).

Nevertheless, Henderson et al. (2005) obtained lower R² values of 0.44 for the surface soil and 0.22 for sub soil prediction of clay content in Australia. The accuracy of prediction depends on the sampling strategy, landscape condition, extent, and resolution. Nevertheless, the validation of these maps show values that are comparable to other studies which used digital soil mapping techniques employing empirical models that relate soil observations and environmental covariates (Liu et al., 2012; McBratney, Mendonça Santos & Minasny, 2003).

To test the accuracy of the predicted water content at −10 kPa (field capacity), we compared the prediction of 7 soil series with an independent water retention dataset. Figure 5A shows that the prediction has an R² = 0.61, and RMSE = 6.3%, which is comparable with values reported by other researchers. Wösten, Pachepsky & Rawls (2001) showed that the RMSE values for the prediction of water content at −10 and −33 kPa based on various studies are between 3 and 7%. Because of water retention is rarely measured in Korean surveys, the validation can only be done on a limited soil samples. Nevertheless, Fig. 5B shows an example of the measured and predicted PAWC for a profile from Songjeong soil series. The profile is a fine loamy, mesic family of Typic Hapludults, with measurement collected from soil core samples from Ap1, B1, and C horizons, located in an apple research orchard field at the College of Agriculture and Life Science, Seoul National University, in Suwon. Although there is variation in the depth of the horizons with the modal profile, the prediction of water content at −10 and −1500 kPa fits reasonably well with the measured data.