Improved runoff simulations for a highly varying soil depth and complex terrain watershed in the Loess Plateau with the Community Land Model version 5

This study aimed to improve runoff simulations and explore deep soil hydrological processes for a watershed in the center of the Loess Plateau (LP), China. This watershed, the Wuding River Basin (WRB), has very complex topography, with soil depths ranging from 0 to 197 m. The hydrological model used for our simulations was Community Land Model version 5 (CLM5) developed by the National Center for Atmospheric Research. Actual soil depths and river channels were incorporated into CLM5 to realistically represent the physical features of the WRB. Through sensitivity tests, CLM5 with 150 soil layers with the observed variable soil depths produced the most reasonable results and was adopted for this study. Our results showed that CLM5 with actual soil depths significantly suppressed unrealistic variations of the simulated subsurface runoff when compared to the default simulations. In addition, when compared with the default version with 20 soil layers, CLM5 with 150 soil layers slightly improved runoff simulations but generated simulations with much smoother vertical water flows that were consistent with the uniform distribution of soil textures in our study watershed. The runoff simulations were further improved by the addition of river channels to CLM5, where the seasonal variability of the simulated runoff was reasonably captured. Moreover, the magnitude of the simulated runoff remarkably decreased with increased soil evaporation by lowering the soil water content threshold, which triggers surface resistance. The lowered threshold was consistent with the loess soil, which has a high sand component. Such soils often generate stronger soil evaporation than soils dominated by clay. Finally, with the above changes in CLM5, the simulated total runoff matched very closely with observations. When compared with those for the default runoff simulations, the correlation coefficient, root mean square error, and Nash–Sutcliffe coefficient for the improved simulations changed dramatically from 0.02, 10.37 mm, and −12.34 to 0.62, 1.8 mm, and 0.61. The results in this study provide strong physical insight for further investigation of hydrological processes in complex terrain with deep soils.

1 Introduction

Understanding runoff processes in regions with very complex topography is important for managing and predicting water resources. Such an understanding can assist in quantifying the allocation of water resources (Chen et al., 2013; Camacho et al., 2015), evaluating surface and groundwater vulnerability to natural and anthropogenic processes (Uhlenbrook et al., 2002), improving drought and flood management (Camacho et al., 2015), and predicting the amount and spatiotemporal distribution of water resources (Saraiva Okello et al., 2018). However, complex topography leads to intricate runoff processes (Jencso et al., 2011), causing uncertain estimation of water resources. Therefore, it is crucial to investigate runoff processes for the wellbeing of topographically complex regions.

As the largest area covered by continuous loess soils in the world (Fu et al., 2017; Zhu et al., 2018), the Loess Plateau (LP) in China has complicated hydrological processes because of its extremely complex topography and unique soil types. Due to an arid and semi-arid climate and a population of more than 100 million (Zhang et al., 2018), this region experiences severe water shortages (Xiao et al., 2019). It is essential to accurately estimate the spatiotemporal distribution of water resources in this region of complex terrain. Soil depth in the LP can reach 350 m (Zhu et al., 2018; Li et al., 2019), making it difficult to measure deep soil hydrological processes and understand runoff generation (Shao et al., 2018; Liu et al., 2012). In addition, terrain in the LP includes loess tablelands, ridges, hills, gullies, and river channels (Fu, 1989), all of which have quite different runoff generation processes (Liu et al., 2012). In loess tablelands with deep water tables (Huang et al., 2013; Shao et al., 2018), the soils store most infiltrated water, generate insignificant surface runoff, and remarkably delay subsurface runoff. Areas in the LP with gullies and river channels usually have high water tables (Liu et al., 2012) and can easily be saturated during precipitation events, generating a large amount of surface runoff. Especially, extreme rainfall events that mostly occur over the summer monsoon season (Tian et al., 2020) produce strong soil erosion and a large amount of fast infiltration-excess surface runoff to the river channels in hillslope areas, sometimes causing severe flooding. In the meantime, the loess soils that dominate the LP have a large capillary porosity, with loose and homogeneous textures due to a high sand component, often resulting in high evaporation (Li et al., 1985; Lei, 1987; Han et al., 1990; Wang and Shao, 2013). A better understanding of the hydrological processes within the complex terrain and special soil types of the LP is vital to improving the prediction of water resources in this region.

Numerical hydrological models are essential tools to investigate runoff processes in the LP. Field measurements such as those from tracer techniques (Huang and Pang, 2011; Li et al., 2017; Huang et al., 2017, 2019; Xiang et al., 2019) have been made to quantify the hydrological processes in the LP, but these measurements have significant limitations, including short temporal and small spatial coverage, which cannot account for the processes at watershed scales. Hydrological models based on mass and energy equations are effective in simulating the long-term spatiotemporal variability of runoff at watershed scales (Döll and Fiedler, 2008; Turkeltaub et al., 2015; Shao et al., 2018). Hydrological models can also simulate the quantity of different components in the water budget (e.g., surface runoff, subsurface, etc.) that are difficult or impossible to be measured directly. Based on detailed soil information at a depth of 98 m at a research site on the Changwu tableland in the LP, Shao et al. (2018) used a hydrological model to generate reasonable simulations for deep soil percolation and groundwater level. Their study provides important clues (e.g., high-resolution soil layering) for exploring deep soil hydrological processes and producing reliable runoff simulations at a watershed scale in the LP. Therefore, it is apparent that hydrological models can overcome the drawbacks of field experiments.

However, in hydrological models, soil depth and river channels are very important in simulating soil water movement and storage and runoff processes, especially in regions with complex topography (Tesfa et al., 2009; Fu et al., 2011). Soil depth is set to a constant in most hydrological models (Shangguan et al., 2017). For example, the Noah (Ek et al., 2003) and Noah-MP (Niu et al., 2011) models have a fixed soil depth of 2 m, which cannot represent the realistic spatial distribution of soil depth in the LP, which ranges from 0 to 350 m. In addition, soil depth in most river channels with exposed bedrock in the LP is close to zero (Jing and Cheng, 1983; Li et al., 2017; Zhu et al., 2018), and areas dominated by these channels are very important in generating runoff. Some hydrological models such as the Community Land Model version 5 (CLM5) and the Soil and Water Assessment Tool (Neitsch et al., 2011) have embedded river routing schemes. In these schemes, the river channels described based on elevation differences still have the same soil depth as other places without these channels, which cannot reflect the actual conditions in the LP and many other regions where soil depth changes significantly across rivers. Thus, soil depth variations and river channels need to be considered in hydrological models for better soil water flow and runoff simulations.

The objective of this study was to use CLM5 to improve runoff simulations and better understand the hydrological processes with varying soil depths for a very complex topography watershed in the LP. To achieve this objective, the highly varying soil depths and river channels were incorporated into CLM5 to realistically represent the features of the watershed. In fact, Brunke et al. (2016) have conducted a study with CLM version 4.5 by including varying soil depths at a global scale where the runoff simulations are focused at grid cell scales, which cannot be evaluated with actual streamflow data. However, evaluating hydrological simulations at watershed scales is essential to improving our understanding of runoff processes. In this study, the most important finding was that river channels where the soil depth is often equal or close to zero played a vital role in runoff simulations especially in complex topography areas. According to our extensive literature search, river channels are not configured in most of existing land surface and hydrological models. In addition, although this study focused on a relatively small watershed, our runoff simulation methods and science ideas can be easily transferred to investigate the hydrological processes in other watersheds across the world with observed soil depth and river channel information. The text is laid out as follows: Sects. 2 and 3 introduce the study area and data, respectively, Sect. 4 provides the model description, Sect. 5 describes the methodology, Sect. 6 includes the results, and the conclusions are in Sect. 7.

2 Study area

The Wuding River Basin (WRB) was selected as the study area. This basin, with an area of about 30 261 km², is in the center of the LP (Fig. 1a), which is the largest continuous loess area in the world (∼ 640 000 km²) (Fu et al., 2017; Zhu et al., 2018). The WRB shows complex geomorphic characteristics including tablelands, ridges, hills, gullies, and river channels (Liu et al., 2012). The main land use types in the WRB are bare ground, grassland, and sparse forest. Across the basin, soil thickness generally ranges from 0 to 200 m (Liu, 2016), and the loess, consisting mainly of silt and sand (Li et al., 1985), is relatively homogeneous in the vertical direction (Huang et al., 2013; Xiang et al., 2019). The WRB has a continental monsoon climate with mean annual precipitation of around 400 mm, about 70 % of which falls during the flood season from June through September, based on observations over the period of 1956–2010 (http://data.cma.cn/, last access: 18 August 2021). Figure 1b shows the geographic distribution of the observed soil depth for the WRB, which is discussed again in Sect. 3.2.

Figure 1 (a) Location of the Loess Plateau (LP) and WRB in China; (b) the geographic distribution of the observed soil depth for the WRB.

3 Data

3.1 Meteorological and runoff data

High-quality meteorological and runoff data for the WRB were used to force and evaluate CLM5, respectively. The Global Soil Wetness Project phase 3 (GSWP3) meteorological dataset (http://hydro.iis.u-tokyo.ac.jp/GSWP3/index.html, last access: 1 July 2020) was selected to drive the model for this study. The GSWP3 dataset contains seven climate forcing variables, including precipitation, air temperature, downward shortwave and longwave radiation, specific humidity, surface pressure, and wind speed. These data cover the period of 1901–2010 with a spatial resolution of 0.5^∘ at a 3 h time step. Meanwhile, we obtained the observed monthly runoff data of the Baijiachuan (BJC) hydrological station from the Data Sharing Network of Earth System Science (http://loess.geodata.cn/index.html, last access: 1 May 2021). The BJC station is located at the WRB outlet and its drainage area covers ∼ 98 % of the basin. These runoff data were used to assess CLM5 output.

3.2 Soil data

Soil depth data for the WRB as shown in Fig. 1b were obtained from different sources. We first collected and recorded 61 soil depths for the WRB and nearby areas from ∼ 15 published papers and books (not cited here). In addition, two soil depth maps for the WRB were obtained from Qi et al. (1991) and Wang (2016) and were digitized. Soil depth data for model grids with gullies and rivers were derived based on digital elevation model (DEM) data. Soil depth in gullies and rivers was assumed to be 0 due to the exposure of bedrock (Jing and Cheng, 1983; Li et al., 2017; Zhu et al., 2018). The elevations of these gully and river channels were retrieved from a DEM at a resolution of 90 m. The differences between these elevations and those at a 5 km resolution were used to represent the soil depth in model grids with gullies and rivers. This is different from river routing that is based only on one DEM. The proportion of the total gully and river area to the entire WRB area (defined as P_gr hereafter) was determined with the Cressman method (Cressman, 1959). A value of 0.3 is suggested by Qi et al. (1991) for the LP. In this study, we identified the optimal P_gr value through sensitivity tests by setting different interpolation radii in the Cressman method. The soil depth data from these sources were then combined and interpolated into a 5 km resolution, still based on the Cressman method.

Soil texture data for the WRB were necessary input into CLM5. These data were derived from a soil type map for the LP (http://loess.geodata.cn, last access: 10 June 2020) and included three soil layers: 0–20, 20–76, and 76–180 cm. For soil layers deeper than 180 cm, the texture data for the 76–180 cm layer were applied.

4 Model description

CLM5 was used in this study for runoff simulations. This model was developed by the National Center for Atmospheric Research. The CLM5 includes one vegetation layer, up to five snow layers, and 20 soil layers. In the model, each grid cell is split into different land units including vegetated surface, lake, urban, glacier, and cropland. The spatial distribution and seasonal climatology of the plant functional types for CLM5 are derived from MODIS satellite land surface data products (Lawrence and Chase, 2007). CLM5 uses the simplified TOPMODEL (Niu et al., 2005) to parameterize runoff, which is partitioned into surface and subsurface runoff. Surface runoff is calculated based on the saturation-excess mechanism. Subsurface runoff is produced when saturated conditions occur within the soil column. CLM5 is attached with a river routing module for runoff simulations. However, in this study, we focused our simulations on a monthly time scale at which the river flow should be able to travel from the farthest point to the outlet of the WRB with an area of 30 261 km² that can easily fit into a 200 by 200 km box. Thus, we turned off the river routing module during our simulations and used the total runoff over the entire watershed for comparison with observations.

In CLM5, soil evaporation is affected by soil resistance, which is associated with a dry surface layer (DSL) (Swenson and Lawrence, 2014). A DSL forms near the soil surface in the model when the soil water content in the top layer is below a threshold value (SWC_th), which is set to 80 % of the soil porosity of the top layer (SWC_sat,1). The formation of the DSL generates soil resistance, limiting soil evaporation. Meanwhile, CLM5 uses Richard's equation and Darcy's law to describe changes in soil water content (SWC) and soil water flux. The soil hydraulic conductivity and retention used in these equations are determined by the soil texture and the SWC of the previous time step, based on Clapp and Hornberger (1978), Cosby et al. (1984), and Lawrence and Slater (2008).

5 Methodology

5.1 Soil layering

As mentioned earlier, actual soil depth in the WRB is strongly variable, with a range of ∼ 0–197 m (Fig. 1b). In our default run, the soil depth in CLM5 was set to a constant of 8.6 m (see Table 2.2.3 in Lawrence et al., 2018) and is discretized into 20 layers defined as hydrological active layers (HALs) to distinguish them from the five bedrock layers set in the model. In this study, we compared the simulations with a default fixed soil depth to those with the observed variable soil depths for the WRB based on the soil depth data shown in Fig. 1b. Eight sensitivity tests were conducted with soil layer numbers (SLNs) of 20, 50, 75, 100, 125, 150, 175, and 200 to determine the optimal soil layering method for runoff simulations in the WRB (Tables 1a and b). In each sensitivity test, the SLN is the same for the entire WRB, and the HAL number is identified based on the input soil depth for each soil column. Layers that are not HALs are treated as bedrock layers and are not used in the hydrology calculations in the model. These sensitivity simulations were compared to those with the default options of CLM to examine how the vertical resolution with observed variable soil depths affected the runoff simulations for the WRB.

Table 1 (a) Thickness (m) of each soil layer for different SLNs (20, 50, 75, and 100). (b) Thickness (m) of each soil layer for different SLNs (125, 150, 175, and 200).

5.2 Model spinup and simulations

All runs in this study needed model spinup to ensure that the soil moisture of each HAL reached equilibrium. We found that the spinup period could last for more than 50 years for different initial SWC conditions and soil depths in the WRB. The initial SWC was set to 0.2 mm³ mm⁻³, and we performed two cycles of continuous simulations over the period of 1901–2010. The first cycle was discarded as spinup, and the second cycle was retained for analysis. Through these spinup runs, the SWC at all model grids can reach the equilibrium state (an example is given in Fig. 5 where the soil has the deepest depth of 197 m in our simulation domain). In this study, each sensitivity run had its own spinup cycles.

Figure 2 Observed monthly precipitation and runoff (black line) and simulated total runoff, surface runoff, and subsurface runoff in the default run from 1956 to 1969. The observed monthly precipitation is for the entire WRB, and the OBS (observations) and simulations are for the BJC hydrological station.

Figure 3 Observed and simulated total runoff, surface runoff, and subsurface runoff from 1956 to 1969 by the run with actual soil depths. The SLN was set to 20.

Figure 4 (a) Observed and simulated monthly runoff for the BJC hydrological station; (b) simulated surface runoff; (c) simulated subsurface runoff; (d) RMSEs of the simulated total runoff. All simulations were produced with different SLN values.

Figure 5 Simulated vertical monthly SWC profiles for the selected point in the WRB during both the spinup (left of the dashed black line) and simulation (right of the dashed black line) periods. All simulations were conducted with different SLN values.

Figure 6 Spatial distributions of river channels (black lines) and soil depths for the WRB with different values of P_gr.

Figure 7 (a) Observed and simulated monthly runoff for the BJC hydrological station; (b) simulated surface runoff; (c) simulated subsurface runoff. All simulations were produced with different P_gr values.

6 Results and analysis

6.1 Default runoff simulation

We conducted a default run to evaluate the performance of the original CLM5 in simulating runoff in the WRB. The model remarkably overestimated monthly total and subsurface runoff when compared with observations from the BJC hydrological station over 1956–1969, a period with minimal human activity (Jiao et al., 2017). The correlation coefficient (R²), root mean square error (RMSE), and Nash–Sutcliffe efficiency (NSE) were 0.02, 10.37 mm, and −12.34, respectively. We can see that the overestimation was due mainly to the unrealistic simulations of subsurface runoff. The reasons for these erroneous simulations are discussed in detail in the following sections.

6.2 Effects of soil depth on runoff simulations

We examined how the simulated runoff for the WRB was affected by the actual soil depths (40–197 m) that were inputted into CLM5 with a default SLN of 20. As shown in Fig. 3, CLM5 with deep soils greatly suppressed the seasonal variability of subsurface runoff and reduced the magnitude of surface runoff when compared with the CLM5 simulations with a uniform soil depth of 8 m. The R², RMSE, and NSE between observations and the simulations with actual soil depths were 0.04, 9.8 mm, and −10.96, respectively. Although the actual soil depth data for the WRB were included in CLM5, the runoff simulations were still remarkably different from observations in both variability and magnitude. Hence, the runoff simulations for the WRB need to be further explored and understood.

6.3 Effects of soil layering on runoff simulations

The eight soil layering methods mentioned in Sect. 3.2 were applied to CLM5 with the actual soil depths for the WRB to investigate the effects of soil layering on the runoff simulations. We can see that all the CLM5 runs generated similar temporal patterns of simulated total runoff, as shown in Fig. 4a. Obviously, the soil layering methods had almost no effect on the surface runoff simulations (Fig. 4b), while these methods did affect the subsurface runoff simulations to some extent (Fig. 4c). When the vertical spatial resolution increased from 20 to 200 soil layers, the RMSE of the simulated total runoff decreased until the SLN was equal to 75, and then the errors reached a minimum for SLN ranging from 100 to 200 (Fig. 4d). Although the model with 75 soil layers seemed to be an efficient case, the soil layering method was further examined with vertical soil moisture profile simulations.

We selected a point (37.53^∘ N, 109.33^∘ E) with the deepest soil depth of 197 m in the WRB to study the soil layering method based on vertical soil moisture profile simulations. As shown in Fig. 5, the coarser-resolution simulations (SLN ≤ 125) resulted in alternating persistent wet–dry layers throughout our study period, and this alternation gradually weakened with increasing SLN. When the SLN was equal to 150, the wet–dry alternation almost disappeared. We examined the model numerical method and found that the coarser resolution numerically caused smaller soil matric potential (SMP) gradients between the soil layers, leading to the wet–dry alternation. These vertical soil moisture simulations indicated that CLM5 could produce smooth soil water flow simulations with at least 150 soil layers at a soil depth of 197 m to avoid these numerical issues, although the RMSE of the simulated total runoff reached the minimum value with SLN equal to 75. Therefore, in the following simulations, we set the model soil layers to 150. With this soil layering, the R², RMSE, and NSE for the total runoff simulations were 0.07, 9.3 mm, and −9.71, respectively.

6.4 Effects of P_gr on runoff simulations

In addition to the actual soil depth and high-resolution soil layering, we prescribed the river channels for the WRB in CLM5 to explore the effects of those channels on runoff simulations. Figure 6 shows the spatial distribution of the river channels for the WRB with different values of P_gr, a proportion of the total river channel area to the entire WRB area, as previously defined. The larger the P_gr, the denser the river channels. Our results showed that CLM5 dramatically improved the simulations of the seasonal variability of total runoff (Fig. 7a), and the R² increased to 0.41–0.56 from 0.07 for the previous simulations. These improvements resulted mainly from the surface runoff simulations with a much higher seasonal variability (Fig. 7b). The subsurface runoff simulations did not show significant changes with the addition of the river channels to CLM5 (Fig. 7c). We can see that CLM5 with P_gr equal to 0.15 produced the lowest RMSE (9.3 mm) and the highest NSE (−9.78), although the R² was not the highest (0.52) with this P_gr value. Moreover, we found that the seasonal peak values of the simulated surface runoff with P_gr values of 0.22 and 0.26 were higher than the observed peak values (figure not shown), which was not realistic. Thus, we selected 0.15 for P_gr for the rest of our simulations.

Table 2 R², RMSE, and NSE for total runoff simulations with different P_gr values.

Figure 8 Time series of observed monthly runoff (black line) for the BJC hydrological station and simulated monthly total (red line), surface (blue line), and subsurface runoff (pink line).

6.5 Water balance analysis

We looked into the water balance for the WRB and attempted to further reduce the biases of the runoff simulations. In the previous sections, the more realistic conditions of the WRB (actual soil depths, high-resolution soil layering, and river channels) were incorporated into CLM5 to improve the runoff simulations, but the simulations were still far away from observations. Tian et al. (2018) indicated that the change in water storage in the WRB approached zero over a period of 13 years. Our study focused on a period of 14 years (1956–1969). Thus, we estimated the mean evapotranspiration (ET) with observed precipitation and runoff over our study period by assuming a water storage change of zero in the WRB as follows:

$$\begin{array}{}\text{(1)}& {\mathrm{ET}}_{\mathrm{avg}}=P_{\mathrm{avg}}-R_{\mathrm{avg}},\end{array}$$

where ET_avg, P_avg, and R_avg are mean ET (mm), precipitation (mm), and runoff (mm) over 1956–1969, respectively. Here, P_avg is 454.7 mm, R_avg is 53.2 mm, and the estimated ET_avg is 401.5 mm. However, the simulated mean ET over the study period was 267.8 mm, which was far below the estimated value. According to the soil evaporation parameterization in CLM5, when the SWC of the top soil layer (SWC₁) was less than SWC_th, a DSL formed to resist soil evaporation. In CLM5, the SWC_th is defined as 80 % of SWC_sat,1. However, previous studies (Lee and Pielke, 1992; Sakaguchi and Zeng, 2009; Flammini et al., 2018) found that soil evaporation starts to decrease significantly when the surface SWC is less than the field capacity. Yang et al. (1985) also found that soil evaporation in the LP slows down when the surface SWC becomes lower than a stable capacity that is close to the field capacity. Thus, in this study, we changed the SWC_th to the SWC_fc,1 to conduct one additional simulation. With this modification, the simulated annual ET fluctuated around the estimated mean ET for our study period (401.5 mm), and the simulated 14-year mean value was 392.5 mm, which was close to the estimated mean. Very importantly, the simulated total runoff drastically reduced to match observations by increasing ET (Fig. 8). When compared with those for the simulations in the last section, R² increased from 0.52 to 0.62, RMSE decreased from 9.3 to 1.8 mm, and NSE increased dramatically from −9.78 to 0.61. Therefore, we remarkably improved runoff simulations with more accurate ET simulations in addition to the more realistic WRB features.

7 Conclusions and discussion

This study was intended to improve runoff simulations with CLM5 for the complex topography of the WRB and to improve our understanding of deep soil hydrological processes. In CLM5, we included actual soil depths for the WRB ranging from 0 to 197 m and added the river channels for this watershed. We tested eight soil layering methods and found that CLM5 with at least 150 soil layers could produce rational simulations for both runoff and the vertical soil moisture profile. Different values of river channel density were examined with CLM5, showing that a ratio of 15 % of the total river channel area to the entire WRB area generated the most reasonable results.

With the above model settings, our simulations showed that CLM5 with actual soil depths greatly suppressed the seasonal variability of simulated subsurface runoff and reduced the simulated surface runoff when compared with the default simulations with a uniform soil depth of 8 m. In addition, CLM5 with finer-resolution soil layering (SLN ≥ 150) led to more accurate runoff and smoother vertical soil water flow simulations than that with coarser-resolution layering, and the latter was consistent with the homogeneous distribution of vertical soil texture in the WRB. The addition of river channels for the WRB to CLM5 significantly increased the seasonal variability of simulated surface runoff, remarkably improving the seasonal variability of simulated total runoff. Moreover, more accurate simulations of soil evaporation in the WRB dramatically reduced the simulated subsurface runoff and improved the total runoff simulations.

Limitations still exist in this study. We used atmospheric forcing data at a 5 km resolution to drive CLM5, but for our study region with very complex terrain, this resolution may not be sufficient and could potentially have generated errors in our simulations. In the meantime, it is very important to expand this study to a larger or even global scale, and accurate soil depth and detailed soil texture data would be vital to such an expanded study. In addition, soil hydraulic properties may change with depth, but this study did not consider such changes, and this needs to be tested in future studies. Despite these limitations, it is clear that our final runoff simulations with an improved CLM5 were highly accurate, and our understanding of deep soil hydrological processes has advanced.