Derivation of the Spatial Distribution of Free Water Storage Capacity Based on Topographic Index

Abstract

Free water storage capacity, an important characteristic of land surface related to runoff process, has a significant influence on runoff generation and separation. It is thus necessary to derive reasonable spatial distribution of free water storage capacity for rainfall-runoff simulation, especially in distributed modeling. In this paper, a topographic index based approach is proposed for the derivation of free water storage capacity spatial distribution. The topographic index, which can be obtained from digital elevation model (DEM), are used to establish a functional relationship with free water storage capacity in the proposed approach. In this case, the spatial variability of free water storage capacity can be directly estimated from the characteristics of watershed topography. This approach was tested at two medium sized watersheds, including Changhua and Chenhe, with the drainage areas of 905 km² and 1395 km², respectively. The results show that locations with larger values of free water storage capacity generally correspond to locations with higher topographic index values, such as riparian region. The estimated spatial distribution of free water storage capacity is also used in a distributed, grid-based Xinanjiang model to simulate 10 flood events for Chenhe Watershed and 17 flood events for Changhua Watershed. Our analysis indicates that the proposed approach based on topographic index can produce reasonable spatial variability of free water storage capacity and is more suitable for flood simulation.

1. Introduction

As hydrological models play increasingly important roles in operational flood forecasting, numerous lumped and physically-based hydrologic models have been developed over the last half-century, including the Standford model in the US, the Tank model in Japan, and the Xinanjiang model in China [1,2,3,4,5,6,7]. Although lumped models have simplistic and effective characteristics, they only partially capture the effect of a watershed’s spatial variability on hydrological processes. Therefore, hydrologists are widely interested in distributed models in order to overcome the deficiencies of lumped models. Distributed models, such as the Grid-Xinanjiang model (GXM), the Systeme Hydrologique Europeen (SHE), Topography based hydrological model (TOPMODEL), and the hydrological enhanced version of the Weather Research and Forecasting model (WRF-Hydro), explicitly resolve spatial variability of topography, land-cover, and rainfall [8,9,10,11,12,13,14]. A majority of distributed hydrologic models are physical-based as they are defined in light of mass and momentum conservation equations. Estimating the spatial distributions of model parameters based on geographical characteristics of watersheds is particularly important for the applications of distributed hydrologic models [15,16,17,18,19,20]. For example, the parameters of MIKE SHE model are mainly dependent on soil properties including soil layer thickness, soil porosity, and hydraulic conductivity [21]. The key parameters of TOPMODEL are related to terrain topography [22].

The GXM, which integrates features from the conceptual rainfall–runoff model and the physically based flow routing model, has been proposed for flood simulation and operational forecasting [23,24,25,26]. The GXM uses identical square cells from the digital elevation model (DEM) as primary computational elements. These elements simulate five hydrologic processes including interception, direct channel precipitation, evapotranspiration, runoff generation, and slope and channel flow routing. One unique feature of the GXM is to parameterize runoff generation due to storage repletion. Runoff is separated into surface, interflow, and groundwater component according to the free water storage. Within each grid cell, runoff is further separated into surface runoff (R_s), interflow (R_i), and groundwater (R_g) components, according to free water storage. For each runoff component, outflow of every grid cell is routed from cell to cell according to the upstream to downstream computation sequencing of grid cells. At the outlet of the watersheds, these runoff components add up to account for the total runoff. The GXM model explicitly parameterizes surface runoff, subsurface flow, and groundwater runoff based on hydraulic routing using the Muskingum method or the diffusion wave method [25]. Keeping the features of the widely used Xinanjiang (XAJ) model facilitates the implementation and operation of the GXM.

Two key parameters of the GXM, as well as the XAJ model, are the spatially variable tension water storage capacity ($W_M$) and free water storage capacity ($S_M$). For W_M, there are physically based methods for estimating its spatial distribution, such as Shi and Yao [25,27]. Yet the free water storage capacity remains empirical. The main goal of this paper, therefore, is to develop a physical based method to estimate the spatial distribution of free water storage capacity.

2. Parameterization of W_M and S_M in the XAJ Model

The XAJ model [28,29,30,31,32] assumes Horton runoff to be the dominant runoff generation mechanism, where runoff, not produced until soil moisture of the aeration zone reaches field capacity (water held by soil against gravity), is equal to rainfall excess without further loss. Since the 1970s, the XAJ model has been successfully applied to numerous humid and semi-humid regions in China. Horton’s runoff has two components: surface and groundwater runoff. Kirkby argued that a third component, interflow, should be considered as infiltration rate varies through soil layer [33,34]. Upper soil often has higher permeability than lower soil leading to a relatively impermeable layer in-between. Infiltrating water from the upper soil will flow laterally, driven by topography—i.e., interflow, which is now included in the XAJ model.

In the XAJ model, a watershed is divided into sub-basins where outflow from each sub-basin is routed down river channels to the outlet of the main watershed. The XAJ model simulations of sub-basin outflow have four major parts: (1) the deficit of soil water storage due to evapotranspiration is distributed through the entire soil layer; (2) runoff production depends on rainfall and soil storage deficit; (3) runoff is partitioned into the surface, subsurface, and groundwater flow; (4) the local runoff is routed to the outlet of a sub-basin

The difference between field capacity and total soil water content is defined as tension water storage capacity, the maximum amount of water available in unsaturated zone [35,36,37,38]. In humid regions, soil moisture often reaches field capacity through the entire soil layer. When soil moisture is below field capacity, infiltration excessive runoff may be produced. Non-uniform spatial distribution of soil water shortage results in a variable distribution of runoff [39,40,41,42,43,44]. A tension water capacity curve is introduced (see Figure 1a) to describe the non-uniform distribution of tension water capacity throughout basin or sub-basin. In Figure 1a, IM is the proportion of impermeable area to total watershed area, AU the initial tension soil water content, P the cumulative precipitation over a certain period of time, K the ratio of potential evapotranspiration to pan evaporation, EM the measured pan evaporation, W the areal mean tension water storage, FR the W-dependent runoff contributing area factor, and B the exponent of the tension water capacity distribution curve. F is the total watershed area. f represents the area of the watershed whose tension water capacity is less than or equal to $W_M$. The tension water capacity at a point, $W_M$, is expressed as,

$$1-\frac{f}{F}={\left(1-\frac{W_M}{W_{M max}}\right)}^B\left(1-IM\right)$$

(1)

where $W_{Mmax}$ is the maximum $W_M$. After soil moisture in an unsaturated zone reaches the field capacity, rainfall entering the reservoir becomes free gravity water until the reservoir is full to produce subsurface runoff. If rain continues, groundwater runoff reaches its limit and the overflow of reservoir becomes surface runoff. The reservoir capacity is assumed to be equal to free water capacity, $S_M$ [45,46,47,48], which plays an important role in the partition of surface and subsurface runoff. Similar to that of tension water capacity, the distribution curve of $S_M$, illustrated in Figure 1b, where the free water capacity, $S_M$, varies from zero to a maximum $S_{Mmax}$:

$$\left(1-\frac{f_r}{FR}\right)={\left(1-\frac{S_M}{S_{M max}}\right)}^{EX}$$

(2)

where $S_{Mmax}$ is the maximum $S_M,$ $f_r$ is the portion of the basin area with free water storage capacity less than or equal to $S_M$, EX the exponent of the free water capacity distribution curve, and FR the runoff contributing area.

3. Spatial Distribution of S_M in the GXM

3.1. Runoff Generation and Separation

Surface runoff $R_s$, interflow $R_i$, and groundwater runoff $R_g$ are related to free water storage S:

$$R_i=K_{i}S$$

(3)

$$R_g=K_{g}S$$

(4)

when $R+S\le S_M$:

$$R_s=0$$

(5)

and when $R+S>S_M$:

$$R_s=R+S-S_M$$

(6)

where $K_i$ is the outflow coefficient of free water storage to interflow, $K_g$ is the outflow coefficient of free water storage to groundwater flow. In the GXM, the sum of K_i and $K_g$ can be fixed and taken as 0.7 according to experience and practical operation in forecasting.

3.2. Estimation of S_M in the GXM

In the GXM, the method of estimating $S_M$, in which the spatial distribution of $S_M$ in the GXM is obtained based on vegetation and soil properties, could be referred to as the “VS method”,

$$S_M=\left({\theta }_s-{\theta }_{fc}\right)L_h$$

(7)

where ${\theta }_s$ is saturated water content, ${\theta }_{fc}$ is field capacity, ${\theta }_{wp}$ is the wilting point, $L_a$ is unsaturated zone thickness, and $L_h$ is humus thickness. $L_a$ and $L_h$ are expressed as,

$$L_a={\epsilon }_{a}Ti+{\epsilon }_b$$

(8)

$$L_h={\vartheta }_k{L}_a$$

(9)

$${\epsilon }_{a}T{i}_{min}+{\epsilon }_b=L_{a max}=\frac{W_{M max}}{{\theta }_{fc,T{i}_{min}}-{\theta }_{wp,T{i}_{min}}}$$

(10)

$${\epsilon }_{a}T{i}_{max}+{\epsilon }_b=L_{a min}=\frac{W_{M min}}{{\theta }_{fc,T{i}_{max}}-{\theta }_{wp,T{i}_{max}}}$$

(11)

where $T{i}_i$ is the topographic index, $T{i}_{min}$ and $T{i}_{max}$ the minimum and maximum topographic index of the watershed, respectively, and $W_{Mmin}$ and $W_{Mmax}$ minimum and maximum $W_M$ of the watershed, respectively. $L_{amin}$ and $L_{amax}$ are the minimum and maximum $L_a$ of the watershed, respectively. The coefficients ${\epsilon }_a$ and ${\epsilon }_b$ are parameterized in terms of topographic index and $W_M$. ${\vartheta }_k$ is the conversion coefficient of humus thickness, depending on unsaturated zone thickness. In practice, the resolutions of vegetation cover and soil type data are often too low for modeling runoff over small watersheds; field observations of ${\theta }_s$, ${\theta }_{fc}$, and ${\theta }_{wp}$ are limited if available at all. Estimation of $S_M$ using high resolution topography data is desirable.

3.3. Estimation of $S_M$ from Topographic Index

The idea of using watershed topography to estimate distributed model parameters has been proposed in previous research, such as TOPMODEL. In TOPMODEL, the terrain index is expressed in terms of terrain topography and hydrological parameters. One parameterization of $W_M$ is the Logarithmic Weibull Curve with Zero Displacement [27].

$$\frac{W_M}{W_{M max}}=EXP\left\{-{\left[\frac{\ln \left(Ti-T{i}_{min}+1\right)}{\alpha }\right]}^{\beta }\right\}$$

(12)

where α and β are empirical coefficients. Combining Equations (1) and (12) leads to a relationship between $S_M$ and topographic index:

$$1-\frac{S_M}{S_{M max}}={\left\{1-EXP{\left[-\left[\frac{\ln \left(Ti-T{i}_{min}+1\right)}{\alpha }\right]\right]}^{\beta }\right\}}^{\raisebox{1ex}{$b$}\!\left/ \!\raisebox{-1ex}{$EX$}\right.}$$

(13)

Equation (13), referred to as “terrain method”, is the proposed method for estimating the spatial distribution of $S_M$, based on topographic characters of watersheds, and establishes a functional relationship between free water storage capacity and topographic index. The $S_{M max}$, B and EX can be roughly estimated according to the application experience of the Xinanjiang model within watersheds. The values of $\alpha$ and $\beta$ can be obtained referring to [27]. $Ti \mathrm{and} T{i}_{min}$ can be calculated according to the geographical characteristics of the watershed.

The Equation (13) show that locations with larger values of free water storage capacity generally correspond to locations with higher topographic index values, such as riparian region. Through this correspondence, the $S_M$ can be estimated by $Ti$ of the known spatial distribution law. For the estimation of S_M spatial distribution, the following steps are taken: (1) calculate topographic index based on DEM; (2) construct Logarithmic Weibull Curve (Equation (12)) with the operate experience of Xinanjiang model and the statistical results of topographic index; (3) based on the Logarithmic Weibull Curve, an equation (Equation (13)) for describing the spatial distribution of S_M is educed by combining Equations (12) and (1). In addition, Equation (13), referred to as “terrain method”, has been implemented in Chenhe and Changhua watersheds to compare the VS method (Equation (7)). The comparison of the XAJ model and GXM with $S_M$, estimated using the terrain method, is discussed in Section 5.

4. Case Studies

4.1. Study Areas and Data

The proposed spatial distribution method was tested on two watersheds. The first was the Chenhe Watershed, located in in the mountainous region of the ShannXi Province. It features a 1395 km² drainage area with nine rain gauges, and elevations ranging from 630 m above sea level at the watershed outlet, to 3747 m (see Figure 2a). The Chenhe Watershed is covered with warm temperate deciduous broad-leaved and conifer forest. The long-term average annual rainfall, pan evaporation, and runoff from 1998 to 2012 were 867 mm, 426 mm, and 471 mm, respectively. Due to the dominance of monsoon climate, more than 60% of annual rainfall occurs during June to September (flood season). The rainfall and discharge records of 10 flood events from the Chenhe Watershed during 2003 to 2012 were used to evaluate the GXM performance, with $S_M$ as a model parameter. The spatial distribution of rainfall was obtained from interpolating the rainfall data from the nine rain gauges (see Figure 2a) using the inverse distance squared procedure. Due to the lack of continuous evaporation data for the calibration period, E-601 pan evaporation data collected at adjacent watershed were used as a surrogate. Other hydro-meteorological data from the Hydrology and Water Resources Survey Bureau of Shannxi Province for the period of 2003–2009 were used for model calibration and validating the remaining data.

The Changhua Watershed is located in the northwest mountainous area of the Zhejiang Province. It features a 905 km² drainage area with elevations ranging from 96 m above sea level at the watershed outlet, to 1748 m (see Figure 2b). With the influence of the summer monsoon, the annual rainfall is as high as 1638 mm. The hourly rainfall and discharge data from 17 flood events (1998–2010) in the Changhua Watershed were obtained from the Hydrology and Water Resources Survey Bureau of the Zhejiang Province. Daily pan evaporation data from the Changhua station are also available.

Two dams of the Huaguangtan Reservoir were built in the Changhua Watershed at 108° E 30° N in 2008. The first dam of 103 m height is the main water retaining structure for water storage and power generation. The second dam, 3 km downstream from the first, is used for water level control related to the operation of the first dam. The water collecting area above the second dam is 350.7 km², which accounts for about 1/3 of the watershed. Hydrological data of the Changhua Watershed are divided into two periods for model calculation: 1998–2008 and 2008–2010. The GXM takes the Huaguangtan Reservoir as a compute node in order to address its effect. In other words, using the Muskingum method, the outflow of the Huaguangtan Reservoir is regarded as the inflow of downstream channel and directly routed to the outlet of watershed. The runoff generation and concentration for the upper region of the reservoir were not taken into consideration for the application of GXM.

The digital elevation model (DEM) with a spatial resolution of 1 km (30″) was obtained from the US Geological Survey [24] and used to derive the topographic attributes of the watershed. Moreover, 1 km of vegetation data are available from the University of Maryland (UMD) Land Cover Classification. Lastly, 10 km of soil texture data are provided by the US Food and Agriculture Organization [24].

4.2. Testing the GXM

The GXM flood simulations used $S_M$ to estimate the terrain method (

Table 1. Statistics of flood simulations.

Period	FloodNo	Robs (mm)	Rsim (mm)	RRE (%)	Qobs (m3/s)	Qsim (m3/s)	RPE (%)	RTE	NS
Cal	2003090319	67.3	70.7	5.07	740	764	3.27	−1	0.96
	2003091711	76.3	79.0	3.54	694	714	2.91	−3	0.89
	2005092413	191.7	223.0	16.30	1191	1295	8.79	−2	0.75
	2005092520	175.1	202.3	15.48	1740	1847	6.15	2	0.78
	2008071909	27.4	29.6	7.90	618	626	1.26	−1	0.79
	2009071419	32.9	38.8	18.04	195	200	2.37	1	0.86
Ver	2010072116	68.9	79.1	14.76	527	624	18.35	−1	0.75
	2011091008	57.4	59.2	2.98	865	805	−6.93	2	0.93
	2011091520	110.1	105.0	−4.56	1200	1131	−5.79	−2	0.96
	2012083013	84.5	78.9	−6.64	1710	1718	0.47	0	0.91

), where FloodNo is the number of flood events, RRE is the relative runoff error, RPE is the relative peak error, RTE the is difference in peak time, and NS is the Nash-Sutcliffe efficiency defined as:

$$RRE=\frac{R_{sim}}{R_{obs}}-1$$

(14)

$$RPE=\frac{ Q_{sim}}{Q_{obs}}-1$$

(15)

$$RTE=T_{sim}-T_{obs}$$

(16)

$$NS=1-\frac{{{\displaystyle \sum }}_{t=1}^z|Q_{sim}^t-Q_{obs}^t{|}^2}{{{\displaystyle \sum }}_{t=1}^z|Q_{sim}^t-\overline{Q_{obs}}{|}^2}$$

(17)

The first four digits of FloodNo represent the year, the fifth and sixth digits represent the month, the seventh and eighth digits represent the day, and the ninth and tenth digits represent the time of flood occurrence. $R_{obs}$ is the observed runoff, $R_{sim}$ the simulated runoff, $Q_{obs}$ the observed peak discharge, $Q_{sim}$ the simulated peak discharge, $T_{sim}$ the simulated time of peak flow, $T_{obs}$ the observed time of peak flow, $Q_{obs}^t$ the observed discharge at time step t, $Q_{sim}^t$ the simulated peak discharge at t, $\overline{Q_{obs}}$ the observed mean over the time period of analysis, and z the total number of times. NS equals to 1 if simulated hydrography agrees with the observed.

5. Results and Discussion

Based on topographic index, the Logarithmic Weibull Curve (Equation (12)) was constructed to derive the equation (Equation (13)) describing the spatial distribution of S_M. Moreover, the application of Equation (13)—referred to as the “terrain method”—was less successful than the XAJ and the GXM in the Chenhe and Changhua watersheds. The spatial distribution of $S_M$ in the Chenhe and Changhua watersheds can be estimated from Ti-$T{i}_{min}$, using the DEM data (Equation (13)). The estimated $S_M$ used in the flood simulations for the Chenhe and Changhua Watershed is shown in Figure 3.

5.1. GXM Simulations of the Chenhe Watershed

The simulation statistics are summarized in Table 1. All calibration and validation simulations were measured using RRE and RPE; the average NS was 0.84 and 0.88, respectively. These GXM simulations justified the topography-based estimates of $S_M$. Three flood events (2003090319, 2003091711, and 2010072116) were selected with the corresponding hydrographs shown in Figure 4. The average free water content changes are shown in Figure 5.

5.2. GXM Simulations for the Changhua Watershed

Six floods were selected over the period of 1998–2008 to calibrate the GXM parameters; four floods used for verification. Over the period of 2008–2010, four floods were selected for calibration, and three floods for validation. The results are shown in

Table 2. Statistics of flood simulations, before the construction of the reservoir.

Period	FloodNo	Robs (mm)	Rsim (mm)	RRE (%)	Qobs (m3/s)	Qsim (m3/s)	RPE (%)	RTE	NS
Cal	1999062404	583.2	530.1	−9.10	2100	1776	−15.43	0	0.91
	1999082923	73.2	72.2	−1.40	950	934	−1.68	1	0.97
	2000053008	61.6	60.3	−2.15	761	780	2.56	0	0.90
	2000060312	44.4	47.9	7.89	548	524	−4.42	0	0.91
	2000062108	78.1	82.5	5.57	700	713	1.89	1	0.86
	2000082407	88.6	95.4	7.73	643	570	−10.48	0	0.90
Ver	2001060923	57.0	63.7	11.71	715	718	0.49	1	0.94
	2002062703	97.4	96.2	−1.22	1340	1382	3.13	1	0.91
	2003051207	72.6	66.7	−8.14	445	456	2.45	1	0.91
	2004051205	87.7	95.2	8.55	368	399	8.37	−2	0.87

and

Table 3. Statistics of flood simulations, after the construction of the reservoir.

Period	FloodNo	Robs (mm)	Rsim (mm)	RRE (%)	Qobs (m3/s)	Qsim (m3/s)	RPE (%)	RTE	NS
Cal	2008061303	46.7	53.5	14.64	507	514	1.47	1	0.90
	2010022508	220.5	234.7	6.41	742	700	−5.72	2	0.89
	2009072316	83.9	100.7	19.93	248	279	12.56	−1	0.93
	2009072805	53.9	52.8	−2.14	499	556	11.47	1	0.84
Ver	2009092917	36.3	32.5	−10.53	580	531	−8.56	−1	0.87
	2009073116	27.7	29.3	5.67	390	350	−10.33	0	0.85
	2010071305	121.3	110.1	−9.30	632	691	9.29	2	0.75

, respectively.

Before the reservoir construction, the model calibrations average NS was 0.91. The model verifications average NS was 0.90. After the reservoir construction, the model calibrations average NS was 0.89. For the model verifications, the average NS is 0.83. The good performance of GXM did not demonstrate that the reservoir construction had no influence on the parameters of runoff generation. Instead, it indicated that the effect of the Huaguangtan Reservoir on the derivation of spatially varied S_M was not obvious under the condition that the reservoir was considered a compute node for the application of GXM.

5.3. Comparison of the XAJ Model and GXM with S_M, Estimated Using the Terrain Method

The GXM flood simulations were compared with the traditional XAJ models in both the Changhua and Chenhe watersheds (see Figure 6). According to the statistics of flood simulations, the XAJ model simulations for the Chenhe and Changhua watersheds before the reservoir construction averaged NS 0.85 and 0.84, respectively. After the reservoir construction in the Changhua Watershed, the average NS becomes 0.81. Both the GXM, with $S_M$ obtained from the terrain index, and the XAJ model, performed well. Furthermore, the GXM performed better for the Changhua Watershed’s 2002062703 flood event, as well as for the Chenhe Watershed’s 2003090319 flood. This enhance performance was due to the non-uniform distribution of rainfall. Figure 7a shows that rainfall during the Chenhe watershe’d 2003090319 flood was mainly concentrated in the upstream area. The water movement of the river took a longer time, causing the delay of the peak flow, when compared to the case of uniform rainfall over the watershed. Figure 6a shows that the peak flow of the XAJ model simulations occured 2 h before the observed; on the contrary, the GXM was in close agreement with the observation. Figure 6b clearly shows that the GXM accurately predicted the peak flow, while the peak flow of the XAJ model simulation was 3 h behind.

5.4. Comparison of Terrain Method and VS Method

In order to simulate the flods in both the Changhua and Chenhe watersheds, we used the VS method to estimate $S_M$ and the terrain method in the GXM. The GXMs average NS—with $S_M$ estimated using the terrain method—was 0.88 and 0.86 for the Changhua and Chenhe watersheds, respectively. The GXMs average NS—with $S_M$ estimated using the VS method—was 0.82 and 0.80 for the Changhua and Chenhe watersheds, respectively. Figure 8 shows that NS was higher when $S_M$ was estimated using the terrain method. When $S_M$ was estimated using the VS method, it was correlated more with soil type than topography, as shown in Figure 9. This caused the spatial discontinuity of $S_M$. Therefore, the terrain method was preferred for the estimation of $S_M$.

The computational efficiency of GXM was of particular importance, given the model was proposed for flood simulation and real-time forecasting. According to the analysis of Yao et al. [49], the cell size of 1 km was more appropriate when the GXM was applied to the medium-sized watersheds (e.g., the Chenhe and Changhua watersheds). In addition, the soil and vegetation data at 1 km resolution were available for our study areas. Therefore, we used the DEM on a grid scale of 1 km. However, the resolution of DEM had an important effect on calculating the topographic index [50,51,52]. The derivation of spatial distribution of S_M was affected by the resolution of the DEM, specifically when using the proposed terrain method. Generally, the homogenization of spatial distribution of topographic index was more obvious when using the coarser DEM. In this case, the heterogeneity of the spatial distribution of S_M derived from the terrain method was smaller. The derived distribution of S_M for the Chenhe and Changhua watersheds may be homogenized, considering that the relatively coarse DEM was used. Furthermore, it will be necessary to investigate whether a more accurate distribution of S_M, as well as an improved performance of GXM, can be obtained through the use of DEM with a higher resolution.

6. Conclusions

In this study, we proposed an improved method to estimate the spatial variability of free water storage capacity, based on topographic index, which was easily extracted from the DEM. The estimated spatial distribution was tested in the GXM. The results indicated that more promising simulations can be obtained from GXM simulations using an original estimation method based on vegetation and soil properties.

We compared the XAJ model and the GXM when using S_M, which was estimated by the terrain method, which simulated the impact of land surface in watersheds after inhomogeneous rainfall during the runoff process. With inhomogeneous rainfall input, the simulation in this study shows that the GXM using S_M—estimated by the terrain method—obtained more accurate results than the XAJ model. Moreover, the comparative analysis of the free water storage capacity of the Changhua and Chenhe watersheds demonstrated how different topographies lead to differences in the spatial distribution characteristics of S_M. Locations with larger values of free water storage capacity generally corresponded to locations with higher topographic index values, such as in the riparian region. In addition, the terrain method overcame the deficiencies of the VS method, as it made S_M continuously change and avoided mutation, due to vegetation or soil changes. This also ensured the continuity of the watersheds hydro-physical characteristics.

Although our preliminary tests demonstrated the potential value of our proposed method, further improvement and research may be needed to achieve a better implementation of this method. Ongoing research could focus on analyzing the effect of different DEM resolutions, such as 30 m and 90 m, on the estimation of the spatial distribution of free water storage capacity, in addition to investigating whether more reasonable spatial patterns of S_M could be obtained using a DEM with a higher resolution.