Abstract
By their nature, wetlands represent an ecosystem base for many concurrent heterogeneous interactions, where the mission of numerical modeling requires a wide range of consistent and reliable datasets from various sources, spatially and temporally. Such a mission usually collides with the existence of tremendous missing in time-series dataset(s). In this context, MIKE SHE was used to construct an integrated surface–subsurface flow model for the Paya Indah wetland in Malaysia where huge gaps exist in the historical datasets of water level and flow rate. To calibrate and validate the model to a satisfactory level, a tri-criteria simulation approach was applied to overcome the occasional missing values in these datasets. This goal was accomplished by calibrating the surface water level and channel flow while simultaneously matching the steady-state subsurface portion of the system wherever water table depth data allowed. Quantitatively, the integrated model scored the highest values of R (0.765 − 0.927) and CE (0.748 − 0.828) during the validation. However, large RMSE values were calculated for the flow rate during calibration at SWL2 (outlet; 0.766) and during validation at Langat River (0.780). This bias was attributed to the low or occasional absence of variation in the historical time-series datasets necessary for the simulation process.
HIGHLIGHTS
The implementation of multi-objective calibration strategy enhances the accuracy of the data-limited wetland model.
Hydro-connectivity of a catchment is the key factor controlling the concept of bridging data gaps.
Plurality/variety of calibration targets can allow data-limited hydrological models to overlap gaps in historical datasets.
This approach provides insights for weighing the hydro-connectivity of a wetland ecosystem.
INTRODUCTION
Depending on the availability, variability and spatial representativeness of the input parameters and on the weight attributed to each hydrologic component of the water cycle, hydrological models as useful tools for water resource management are meant to demonstrate satisfactorily accurate simulation of some or all of the hydrological processes occurring within a watershed (Nguyen & Dietrich 2018; Parra et al. 2019; Bizhanimanzar et al. 2020). In fact, wetland ecosystems worldwide are generally not understood adequately owing to their complex ecohydrological nature, whereby periodic fluctuations of surface water and groundwater head increase the complexity of the hydrological system of such watersheds (Pechlivanidis et al. 2011; Wu et al. 2014a, 2014b; House et al. 2016; Beaulieu et al. 2018; Mehdi et al. 2018). The functionality of using multiple criteria instead of single calibration target for hydrological models, which was established using different modeling codes and a variety of input parameters, depends on the quality and representativeness of the calibration datasets (Pechlivanidis et al. 2011; Wu & Liu 2012; Liu et al. 2016; Nguyen & Dietrich 2018; Bizhanimanzar et al. 2020). In most cases, such models failed to represent precisely the dynamic water level in shallow aquifers owing to unsatisfactory base-flow simulation and oversimplified parameterization of the unsaturated flow which eventually led to replacement of conceptual groundwater modules with external coupled surface water–groundwater models (Alminagorta et al. 2016; Mehdi et al. 2018; Nguyen & Dietrich 2018; Adeyeri et al. 2020; Bizhanimanzar et al. 2020).
The multi-criteria calibration approach was found to have satisfactory performance for catchment hydrological models. Such models, which include MIKE SHE and SWAT, require a multi-objective calibration process to compensate scarcity in both temporal and spatial datasets and to overcome model limitations that might produce biased output (Liu et al. 2016; Ji et al. 2019). Even though adjustments between several criteria are crucial to balance conflicting objectives, it remains difficult to identify an optimal dataset; consequently, the multi-criteria simulation approach is required to consider different calibration targets concurrently (Pechlivanidis et al. 2011; Ahmadi et al. 2014; Wu & Liu 2014; Wu et al. 2014a, 2014b; Ji et al. 2019).
The full-integration criterion requires a wide range of input parameters, high heterogenetic levels of spatial distribution and a linear time-varying system. This, in turn, might dramatically increase the chance of occurrence of some major calibration issues such as over-parameterization, equifinality and model structure uncertainty, the very thing that raises the threshold of acceptability for users (Ma et al. 2016). Many attempts have been made to overcome the limitations of MIKE SHE model performance, including addressing calibration issues through the inclusion of an uncertainty analysis method that enabled the model to solve and satisfactorily interpret the entire range of ecohydrological processes (Refsgaard et al. 2010, 2012; Ewen et al. 2012).
Nevertheless, the major challenges linked to MIKE SHE model calibration, especially for medium- to large-scale catchments, might derive from the complex nature of the study site in question and limitations of the input data (Vansteenkiste et al. 2013; Choi et al. 2015; Sandu & Virsta 2015; Ma et al. 2016; Prucha et al. 2016). In such situations, for optimal model performance, it is crucial to curb unreliable judgments by maintaining similarities between the model and reality during the calibration processes using different sets of targeted parameters (Choi et al. 2015; Sandu & Virsta 2015; Paparrizos & Maris 2017).
By their nature, wetlands are a complex ecohydrological system where the mission of numerical modeling requires a wide range of consistent and reliable time-series datasets from a variety of different sources. However, such a mission usually collides with the existence of huge gaps in historical datasets, which in turn undermine the key processes of model performance evaluation, i.e., calibration and validation. In practice, routine collection of hydrometeorological data is rare in many wetlands or watersheds because of the problem of vandalism. Consequently, measurements in such regions are almost invariably limited to a single site or, at most, to a small number of locations. Accordingly, the spatial variability of hydrological conditions might not be maintained, which is the thing that most affects the satisfactory application of distributed hydrological models to the simulation of spatially variable processes. Furthermore, in wetland hydrological models, data limitations mostly necessitate (1) the application of a simple approach, e.g., interpolation for rainfall and observed time-series water level fluctuations or (2) the adoption of certain uniform values, e.g., for potential evapotranspiration. However, these options can result in overall bias in model predictivity, especially when huge discontinuity is encountered in the observational dataset of one or many of the calibration targets. The approach of adopting multiple calibration targets allows better use of the available input dataset, regardless of its consistency, for simulation of hydrological processes by maintaining a satisfactory level of conformity between the different targets through the statistical measures embedded in the model. Thus, the auto-calibration process becomes most effective when there is variability of the calibration targets. This study, therefore, aimed to overcome limitations in observational data continuity by increasing the number of calibration targets, as well as the number of simulation points for each target to be validated later using accurate field datasets.
Study area
Input datasets
Station . | Weighted factor (%100) . |
---|---|
Sungai Manggis | 40.10 |
Bukit Cheeding | 36.17 |
Perang Besar | 20.78 |
KLIA | 2.95 |
Station . | Weighted factor (%100) . |
---|---|
Sungai Manggis | 40.10 |
Bukit Cheeding | 36.17 |
Perang Besar | 20.78 |
KLIA | 2.95 |
Potential evapotranspiration (ETp) was estimated as a function of ET0 by means of a crop coefficient, i.e., ETp = ET0 × kc, where kc is the crop coefficient (FAO 2012). The MIKE SHE code estimates the actual evapotranspiration (ETa) based on the model of Kristensen & Jensen (1975). The MIKE SHE model required the leaf area index (LAI) and the root mass distribution (RD) for each land use type. The vegetation-specific parameters for the PIW were retrieved from the land use map of Kuala Langat for the years 1998 and 2006. For the PIW catchment, the total number of land used types identified was 10 (Figure 2(d); Table 2). The distribution of vegetation parameters was based on identification of the dominant characteristic vegetation and land use types within the basin. The main impact on land use has been logging activities, especially in the area to the north of Main Lake where the forest has been removed and the land left barren. Such activity has been the primary cause of the decay and subsidence of the peatland blanket in this area.
Vegetation development . | Township . | Citrus . | ||||||
---|---|---|---|---|---|---|---|---|
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 2 | 200 | 0.50 | 0 | 4.5 | 1,250 | 1 |
Harvest | 365 | 2 | 200 | 0.50 | 365 | 4.5 | 1,250 | 1 |
Vegetation development . | Sugarcane . | Pasture . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0.0 | 1.0 | 500 | 1 | 0 | 3 | 600 | 1 |
Early season | 60 | 2.0 | 1,000 | 1 | 150 | 3 | 600 | 1 |
Mid season | 90 | 2.5 | 1,500 | 1 | 240 | 4 | 600 | 1 |
Late season | 120 | 3.5 | 1,500 | 1 | 330 | 4 | 600 | 1 |
1st harvest | 150 | 4.5 | 1,500 | 1 | 365 | 3 | 600 | 1 |
2nd harvest | 210 | 5.5 | 1,500 | 1 | ||||
3rd harvest | 365 | 6.0 | 1,500 | 1 | ||||
Vegetation development . | Truck crop . | Grass . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 4.5 | 750 | 1 | 0 | 3 | 750 | 1 |
Harvest | 365 | 4.5 | 750 | 1 | 365 | 3 | 750 | 1 |
Vegetation development . | Shrub . | Marsh . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 3 | 750 | 1 | 0 | 2 | 750 | 1 |
Harvest | 365 | 3 | 750 | 1 | 150 | 3 | 750 | 1 |
365 | 3 | 750 | 1 | |||||
Vegetation development . | Sparse forest . | Oil palm . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 3 | 750 | 1 | 0 | 2.5 | 1 | 0.85 |
Harvest | 365 | 3 | 750 | 1 | 365 | 2.5 | 1 | 0.90 |
Vegetation development . | Township . | Citrus . | ||||||
---|---|---|---|---|---|---|---|---|
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 2 | 200 | 0.50 | 0 | 4.5 | 1,250 | 1 |
Harvest | 365 | 2 | 200 | 0.50 | 365 | 4.5 | 1,250 | 1 |
Vegetation development . | Sugarcane . | Pasture . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0.0 | 1.0 | 500 | 1 | 0 | 3 | 600 | 1 |
Early season | 60 | 2.0 | 1,000 | 1 | 150 | 3 | 600 | 1 |
Mid season | 90 | 2.5 | 1,500 | 1 | 240 | 4 | 600 | 1 |
Late season | 120 | 3.5 | 1,500 | 1 | 330 | 4 | 600 | 1 |
1st harvest | 150 | 4.5 | 1,500 | 1 | 365 | 3 | 600 | 1 |
2nd harvest | 210 | 5.5 | 1,500 | 1 | ||||
3rd harvest | 365 | 6.0 | 1,500 | 1 | ||||
Vegetation development . | Truck crop . | Grass . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 4.5 | 750 | 1 | 0 | 3 | 750 | 1 |
Harvest | 365 | 4.5 | 750 | 1 | 365 | 3 | 750 | 1 |
Vegetation development . | Shrub . | Marsh . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 3 | 750 | 1 | 0 | 2 | 750 | 1 |
Harvest | 365 | 3 | 750 | 1 | 150 | 3 | 750 | 1 |
365 | 3 | 750 | 1 | |||||
Vegetation development . | Sparse forest . | Oil palm . | ||||||
D.O.G . | LAI . | RD . | Kc . | D.O.G . | LAI . | RD . | Kc . | |
Planting | 0 | 3 | 750 | 1 | 0 | 2.5 | 1 | 0.85 |
Harvest | 365 | 3 | 750 | 1 | 365 | 2.5 | 1 | 0.90 |
LAI, leaf area index (m2 m−2); RD, rooting depth (mm); Kc, crop coefficient (–); D.O.G, days of growth.
Hydrological process
The MIKE SHE model describes hydrological processes mostly based on the laws of conservation of mass, momentum and energy. The 1D and 2D diffusive wave Saint Venant equations are used to describe channel and overland flow, respectively. The methods of Kristensen and Jensen (Refsgaard et al. 2010) are used for evapotranspiration, the 1D equation of Richards (Richards 1933) is used for the unsaturated zone flow and a 3D Boussinesq equation is used for the saturated zone flow (Boussinesq 1904). These partial differential equations are solved using finite difference methods, while other methods (namely, interception/evapotranspiration and snowmelt) in the model are based on empirical equations obtained from independent experimental research (DHI 2004).
Model domain and discretization
The area was discretized into a mesh of 200 × 200 m comprising 6,500 cells. The total number of computational cells in the groundwater model was approximately 19,500. The hydraulic model is overlain by both the unsaturated zone model and the overland flow model. The model uses Cassini coordinates of Selangor State in Malaysia in metric units. Despite the relatively coarse discretization of 200 × 200 m, it appears suitable for calculating both the overall water balance and the impact of large-scale developments (Vázquez & Feyen 2007; Im et al. 2009; Refsgaard et al. 2012) such as in the town of Cyberjaya adjacent to the northeast corner of the catchment (Figure 2(a)). However, the model could easily be refined or zoomed-in models created.
Overland flow and detention storage
Overland sheet flow occurs when the water depth on the ground surface exceeds 2 mm. Other than surface topographic elevation, detention storage and the roughness coefficient (Manning's number; N), the required input data for the surface water model (MIKE11 model) consists of branch networks, cross sections, control structure geometry and operation schedules. The model area is dominated by reasonably flat peat swamp with little topographic relief. Therefore, overland flow was considered of little importance in the peat swamp where the undulating ground surface restricts surface flow. However, for the other types of soils, the overland flow direction and velocity were determined based on the ground surface slope and the roughness coefficient. The dimensionless values of Manning's number assigned for the flow rate in PIW catchment ranged between ‘5’ for the peat swamp and ‘10’ for channels.
Flooded areas
To describe flow attenuation and surface water storage, information on floodplain dynamics was crucial for simulation of river/floodplain interaction. Accordingly, with the exception of the narrow tongue at the northeastern corner of the PIW catchment, the ground surface of the Kuala Langat swamp forest was generally considered flat with a number of flooded areas (lakes), floodplains and depressions. The Paya Indah lake system was assigned as a flooded area with a code given for each lake (Figure 2(e)), whereas the depression and floodplain data were derived from the topographical model by comparing riverbank elevations with surrounding surface elevations.
Hydrogeological parameters
Hydrogeologically, the PIW model comprises three geological layers which are, from top to bottom, peat (≈8 m), silty-clayey sand (≈20 m) and silty-sandy gravel (≈30 m). Hydraulic properties for each layer including hydraulic conductivity, transmissivity and storage and leakage coefficients datasets were used to develop the geological model (Figure 2(f)).
The low hydraulic conductivity of Layer 2 (aquitard) indicates that the PIW lake system is not in direct hydraulic contact with the aquifer, which in turn proves that the lakes do not recharge the aquifer. This conclusion is supported by the fact that these lakes have been filled with mining material (clays/silts) that is characterized by very low hydraulic conductivity (9.95 × 10−7 to 1.2 × 10−8 m/s). Furthermore, the large hydraulic gradients observed in the Paya Indah lake system (different lake water levels) would not be possible under natural conditions if the lakes were settled in highly permeable soils, especially during periods of little or no inflow to Paya Indah. The horizontal conductivity of the aquifer ranges between 4.1 × 10−3 and 3.3 × 10−5 m/s, whereas the vertical conductivity estimated to be (0.1) × (value of horizontal conductivity measured for each layer).
The storage coefficient ranges from 0.001 to 4.5 × 10−5. Withdrawal of groundwater was specified at the Megasteel well-field. For base-flow estimation, the value of 0.001 was assigned for the specific yield of the shallow aquifer of Layer 1 (unconfined aquifer). For Layer 3 (deep aquifer), a storage coefficient value of 4.5 × 10−5 was assigned.
METHODS
Water movement in the unsaturated zone
Assessment of calibration procedure
Given the complexity and data-limited situation of the PIW catchment reflected in the calibration procedure of the distributed parameters, the calibration process aimed to obtain a set of model parameters that could provide satisfactory agreement between model results and field observations, i.e., the calibration process takes more or fewer iterations depending on how well the model performs. In fact, automatic methods in MIKE SHE allow for quick and objective model calibration. However, these methods did not provide satisfactory results obviously due to big gaps in the calibration target datasets. We found that combining automatic with some manual adjustment for the sensitive parameters of PIW is the best approach to handle the complexity of the PIW model. This approach provides a stronger check on model performance than a single auto-calibration method that accumulated model errors over a large range of hydrologic behaviors in the PIW model. In this context, three main groups of input datasets were selected to control the accuracy of the calibration process by adjusting the selected parameter according to its specific measured/selected range. The best values and ranges were obtained through the highest performing model (Table 3), which was evaluated quantitatively through a set of statistical measures embedded in the MIKE SHE system (Tables 4 and 5). The model performance was also assessed visually to determine the degree of conformity between the input data and the modeled values. These inputs included the optimal values of (1) the hydraulic conductivity (Ks) in vadose zone and evapotranspiration rates that were used in the two-layer water balance of the unsaturated zone model to ensure correct simulation of the water yield volume, (2) Manning's roughness coefficient for overland grid cells and channel grid cells and (3) the saturated hydraulic conductivity, i.e., the specific yield of the aquifer of the groundwater model and the model for aquifer–channel water exchange (Table 3). Depending on the availability of the input time-series datasets required for a successful calibration process, the model was calibrated against the surface water level, groundwater head and channel flow rate for the period from 1 July 1999 to 31 October 2004.
Module . | Sensitive parameter . | Optimal values and ranges . | |
---|---|---|---|
Evapotranspiration | Interception coefficient Cint | 0.5 mm | |
C1 | 0.3 | ||
C2 | 0.2 | ||
C3 | 30 mm day−1 | ||
Aroot | 1,250 mm | ||
1,000 mm | |||
700 mm | |||
600 mm | |||
Overland flow (peat swamp) | Manning's number (n) m−1/3 s−1 | 5 | |
Channel flow | Manning's number (n) m−1/3 s−1 | 10 | |
Unsaturated zone: Soil series of Layer 2 (silty clay formation) | Hydraulic conductivity (Ks) | Mined and soil association | 4.42 × 10−8 (m/s) |
Serdang-Bungor Munchong | 4.53 × 10−7 (m/s) | ||
Selangor-Kanchung | 3.95 × 10−7 (m/s) | ||
Perang | 6.63 × 10−7 (m/s) | ||
Saturated zone | Peat aquifer (Layer 1): (clayey peat) | Horizontal hydraulic conductivity (Ks) | 1.0 × 10−4–7.118 × 10−5 (m/s) |
Vertical hydraulic conductivity (Ks) | 1.0 × 10−5–7.12 × 10−6(m/s) | ||
Deep aquifer (Layer 3): (silty sand) | Horizontal hydraulic conductivity (Ks) | 1.0 × 10−5–7.0 × 10−5 (m/s) | |
Vertical hydraulic conductivity (Ks) | 1.0 × 10−5–7.0 × 10−6 (m/s) |
Module . | Sensitive parameter . | Optimal values and ranges . | |
---|---|---|---|
Evapotranspiration | Interception coefficient Cint | 0.5 mm | |
C1 | 0.3 | ||
C2 | 0.2 | ||
C3 | 30 mm day−1 | ||
Aroot | 1,250 mm | ||
1,000 mm | |||
700 mm | |||
600 mm | |||
Overland flow (peat swamp) | Manning's number (n) m−1/3 s−1 | 5 | |
Channel flow | Manning's number (n) m−1/3 s−1 | 10 | |
Unsaturated zone: Soil series of Layer 2 (silty clay formation) | Hydraulic conductivity (Ks) | Mined and soil association | 4.42 × 10−8 (m/s) |
Serdang-Bungor Munchong | 4.53 × 10−7 (m/s) | ||
Selangor-Kanchung | 3.95 × 10−7 (m/s) | ||
Perang | 6.63 × 10−7 (m/s) | ||
Saturated zone | Peat aquifer (Layer 1): (clayey peat) | Horizontal hydraulic conductivity (Ks) | 1.0 × 10−4–7.118 × 10−5 (m/s) |
Vertical hydraulic conductivity (Ks) | 1.0 × 10−5–7.12 × 10−6(m/s) | ||
Deep aquifer (Layer 3): (silty sand) | Horizontal hydraulic conductivity (Ks) | 1.0 × 10−5–7.0 × 10−5 (m/s) | |
Vertical hydraulic conductivity (Ks) | 1.0 × 10−5–7.0 × 10−6 (m/s) |
Elements of multi-objective calibration strategy . | Coordinates (Cassini system) . | Evaluation criteriaa . | PIW prediction variablesb . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
Calibration target . | Category . | Name of calibration point . | X . | Y . | MAE . | RMSEc . | STDres . | R . | CE . | Slope . | Intercept . | r2 . |
Water level (m) | Surface water | Inlet (SWL1) | −7,770.00 | −33,950.0 | 0.167 | 0.217 | 0.207 | 0.848 | 0.812 | |||
Visitor Lake | −7,373.50 | −34,496.3 | 0.063 | 0.076 | 0.076 | 0.753 | 0.737 | |||||
Main Lake | −9,170.00 | −33,280.0 | 0.041 | 0.052 | 0.052 | 0.755 | 0.740 | |||||
Driftwood Lake | −9,562.21 | −33,496.5 | 0.084 | 0.097 | 0.097 | 0.960 | 0.846 | |||||
Perch Lake | −9,958.52 | −34,811.5 | 0.071 | 0.086 | 0.085 | 0.757 | 0.739 | |||||
Marsh Lake | −9,733.34 | −34,415.2 | 0.079 | 0.084 | 0.066 | 0.806 | 0.778 | |||||
Crocodile Lake | −8,620.00 | −34,720.0 | 0.060 | 0.075 | 0.075 | 0.776 | 0.726 | |||||
Hippo Lake | −8,620.00 | −34,720.0 | 0.103 | 0.126 | 0.107 | 0.764 | 0.751 | |||||
Chalet Lake | −8,280.00 | −35,300.0 | 0.071 | 0.098 | 0.097 | 0.785 | 0.764 | |||||
Typha Lake (Lotus T)d | −7,346.48 | −35,460.1 | 0.115 | 0.143 | 0.142 | 0.653 | 0.613 | N/A | N/A | N/A | ||
Lotus Lake (Lotus L) | −7,346.48 | −35,460.1 | 0.124 | 0.149 | 0.142 | 0.645 | 0.598 | |||||
Lotus-Outlet (SWL2) | −11,690.0 | −36,140.0 | 0.170 | 0.205 | 0.205 | 0.458 | 0.425 | |||||
Groundwater | BH1 | −10,094.3 | −26,813.7 | 0.144 | 0.180 | 0.172 | 0.824 | 0.782 | ||||
BH2 | −8,764.66 | −31,583.4 | 0.134 | 0.175 | 0.161 | 0.793 | 0.761 | |||||
BH3 | −9,253.24 | −32,662.5 | 0.073 | 0.091 | 0.091 | 0.866 | 0.791 | |||||
BH4 | −7,592.13 | −33,242.2 | 0.082 | 0.105 | 0.096 | 0.752 | 0.733 | |||||
BH5 | −11,028.3 | −37,772.6 | 0.199 | 0.236 | 0.213 | 0.824 | 0.758 | |||||
BH6 | −13,186.1 | −41,540.7 | 0.253 | 0.317 | 0.316 | 0.693 | 0.634 | |||||
BH7 | −16,698.5 | −39,235.0 | 0.181 | 0.222 | 0.170 | 0.624 | 0.587 | |||||
BH8 | −17,616.6 | −38,289.4 | 0.095 | 0.117 | 0.116 | 0.659 | 0.603 | |||||
Channel flow (m3/s) | Surface water | Inlet | −7,770.00 | −33,950.0 | 0.005 | 0.090 | 0.086 | 0.880 | 0.808 | 0.93** | 0.16** | 0.78 |
Lotus-Outlet | −11,690.0 | −36,140.0 | −0.094 | 0.766 | 0.743 | 0.570 | 0.430 | 0.35 | 0.27* | 0.21 |
Elements of multi-objective calibration strategy . | Coordinates (Cassini system) . | Evaluation criteriaa . | PIW prediction variablesb . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
Calibration target . | Category . | Name of calibration point . | X . | Y . | MAE . | RMSEc . | STDres . | R . | CE . | Slope . | Intercept . | r2 . |
Water level (m) | Surface water | Inlet (SWL1) | −7,770.00 | −33,950.0 | 0.167 | 0.217 | 0.207 | 0.848 | 0.812 | |||
Visitor Lake | −7,373.50 | −34,496.3 | 0.063 | 0.076 | 0.076 | 0.753 | 0.737 | |||||
Main Lake | −9,170.00 | −33,280.0 | 0.041 | 0.052 | 0.052 | 0.755 | 0.740 | |||||
Driftwood Lake | −9,562.21 | −33,496.5 | 0.084 | 0.097 | 0.097 | 0.960 | 0.846 | |||||
Perch Lake | −9,958.52 | −34,811.5 | 0.071 | 0.086 | 0.085 | 0.757 | 0.739 | |||||
Marsh Lake | −9,733.34 | −34,415.2 | 0.079 | 0.084 | 0.066 | 0.806 | 0.778 | |||||
Crocodile Lake | −8,620.00 | −34,720.0 | 0.060 | 0.075 | 0.075 | 0.776 | 0.726 | |||||
Hippo Lake | −8,620.00 | −34,720.0 | 0.103 | 0.126 | 0.107 | 0.764 | 0.751 | |||||
Chalet Lake | −8,280.00 | −35,300.0 | 0.071 | 0.098 | 0.097 | 0.785 | 0.764 | |||||
Typha Lake (Lotus T)d | −7,346.48 | −35,460.1 | 0.115 | 0.143 | 0.142 | 0.653 | 0.613 | N/A | N/A | N/A | ||
Lotus Lake (Lotus L) | −7,346.48 | −35,460.1 | 0.124 | 0.149 | 0.142 | 0.645 | 0.598 | |||||
Lotus-Outlet (SWL2) | −11,690.0 | −36,140.0 | 0.170 | 0.205 | 0.205 | 0.458 | 0.425 | |||||
Groundwater | BH1 | −10,094.3 | −26,813.7 | 0.144 | 0.180 | 0.172 | 0.824 | 0.782 | ||||
BH2 | −8,764.66 | −31,583.4 | 0.134 | 0.175 | 0.161 | 0.793 | 0.761 | |||||
BH3 | −9,253.24 | −32,662.5 | 0.073 | 0.091 | 0.091 | 0.866 | 0.791 | |||||
BH4 | −7,592.13 | −33,242.2 | 0.082 | 0.105 | 0.096 | 0.752 | 0.733 | |||||
BH5 | −11,028.3 | −37,772.6 | 0.199 | 0.236 | 0.213 | 0.824 | 0.758 | |||||
BH6 | −13,186.1 | −41,540.7 | 0.253 | 0.317 | 0.316 | 0.693 | 0.634 | |||||
BH7 | −16,698.5 | −39,235.0 | 0.181 | 0.222 | 0.170 | 0.624 | 0.587 | |||||
BH8 | −17,616.6 | −38,289.4 | 0.095 | 0.117 | 0.116 | 0.659 | 0.603 | |||||
Channel flow (m3/s) | Surface water | Inlet | −7,770.00 | −33,950.0 | 0.005 | 0.090 | 0.086 | 0.880 | 0.808 | 0.93** | 0.16** | 0.78 |
Lotus-Outlet | −11,690.0 | −36,140.0 | −0.094 | 0.766 | 0.743 | 0.570 | 0.430 | 0.35 | 0.27* | 0.21 |
aEvaluation criteria: MAE, mean absolute error; RMSE, root mean square error; STDres, standard deviation of the residuals; R, correlation.
CE: Nash and Sutcliffe coefficient of efficiency.
bConsistency of overall flow of PIW model predictions beyond its development conditions.
cRMSE units are (m) for surface water and groundwater levels and (m3/s) for channel flow rate.
dTypha Lake is considered an extension of Lotus Lake.
N/A, not applicable.
*Significant slope values (P ≤ 0.05).
**Very significant slope values (P ≤ 0.01).
Elements of multi-objective validation strategy . | Coordinates (Cassini system) . | Evaluation criteriaa . | PIW prediction variablesb . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
Validation target . | Category . | Name of validation point . | X . | Y . | MAE . | RMSEc . | STDres . | R . | CE . | Slope . | Intercept . | r2 . |
Water level (m) | Surface water | Visitor Lake | −7,373.50 | −34,496.3 | 0.077 | 0.091 | 0.083 | 0.771 | 0.763 | |||
Main Lake | −9,170.00 | −33,280.0 | 0.056 | 0.073 | 0.069 | 0.765 | 0.780 | |||||
Tin Lake | −9,733.34 | −34,415.2 | 0.049 | 0.057 | 0.053 | 0.814 | 0.794 | |||||
Crocodile Lake | −8,620.00 | −34,720.0 | 0.063 | 0.076 | 0.075 | 0.792 | 0.751 | |||||
Hippo Lake | −8,620.00 | −34,720.0 | 0.037 | 0.049 | 0.048 | 0.783 | 0.776 | N/A | N/A | N/A | ||
Chalet Lake | −8,28000.00 | −35,300.0 | 0.091 | 0.105 | 0.101 | 0.807 | 0.782 | |||||
Typha Lake (Lotus T )d | −7,346.48 | −35,460.1 | 0.070 | 0.087 | 0.078 | 0.871 | 0.760 | |||||
Lotus Lake (Lotus L) | −7,346.48 | −35,460.1 | 0.068 | 0.088 | 0.081 | 0.866 | 0.748 | |||||
Groundwater | BH3 | −9,253.24 | −32,662.5 | 0.101 | 0.132 | 0.132 | 0.791 | 0.774 | ||||
BH5 | −11,028.3 | −37,772.6 | 0.058 | 0.071 | 0.053 | 0.683 | 0.667 | |||||
Channel flow (m3/s) | Surface water | Langat River | −2,207.72 | −35,100.4 | 0.480 | 0.780 | 0.770 | 0.877 | 0.772 | 0.95** | 0.16** | 0.79 |
Lotus-Outlet | −11,690.0 | −36,140.0 | 0.031 | 0.017 | 0.016 | 0.927 | 0.828 | 0.90** | 0.07** | 0.86 |
Elements of multi-objective validation strategy . | Coordinates (Cassini system) . | Evaluation criteriaa . | PIW prediction variablesb . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
Validation target . | Category . | Name of validation point . | X . | Y . | MAE . | RMSEc . | STDres . | R . | CE . | Slope . | Intercept . | r2 . |
Water level (m) | Surface water | Visitor Lake | −7,373.50 | −34,496.3 | 0.077 | 0.091 | 0.083 | 0.771 | 0.763 | |||
Main Lake | −9,170.00 | −33,280.0 | 0.056 | 0.073 | 0.069 | 0.765 | 0.780 | |||||
Tin Lake | −9,733.34 | −34,415.2 | 0.049 | 0.057 | 0.053 | 0.814 | 0.794 | |||||
Crocodile Lake | −8,620.00 | −34,720.0 | 0.063 | 0.076 | 0.075 | 0.792 | 0.751 | |||||
Hippo Lake | −8,620.00 | −34,720.0 | 0.037 | 0.049 | 0.048 | 0.783 | 0.776 | N/A | N/A | N/A | ||
Chalet Lake | −8,28000.00 | −35,300.0 | 0.091 | 0.105 | 0.101 | 0.807 | 0.782 | |||||
Typha Lake (Lotus T )d | −7,346.48 | −35,460.1 | 0.070 | 0.087 | 0.078 | 0.871 | 0.760 | |||||
Lotus Lake (Lotus L) | −7,346.48 | −35,460.1 | 0.068 | 0.088 | 0.081 | 0.866 | 0.748 | |||||
Groundwater | BH3 | −9,253.24 | −32,662.5 | 0.101 | 0.132 | 0.132 | 0.791 | 0.774 | ||||
BH5 | −11,028.3 | −37,772.6 | 0.058 | 0.071 | 0.053 | 0.683 | 0.667 | |||||
Channel flow (m3/s) | Surface water | Langat River | −2,207.72 | −35,100.4 | 0.480 | 0.780 | 0.770 | 0.877 | 0.772 | 0.95** | 0.16** | 0.79 |
Lotus-Outlet | −11,690.0 | −36,140.0 | 0.031 | 0.017 | 0.016 | 0.927 | 0.828 | 0.90** | 0.07** | 0.86 |
aEvaluation criteria: MAE, mean absolute error; RMSE, root mean square error; STDres, standard deviation of the residuals; R, correlation.
CE: Nash and Sutcliffe coefficient of efficiency.
bConsistency of overall flow of PIW model predictions beyond its development conditions.
cRMSE units are (m) for surface water and groundwater levels and (m3/s) for channel flow rate.
dTypha Lake is considered an extension of Lotus Lake.
N/A, not applicable.
*Significant slope values (P ≤ 0.05).
**Very significant slope values (P ≤ 0.01).
Assessment of validation procedure
The validation process covered the period from 1 August 2007 to 1 August 2008. The surface water level and groundwater head were validated at only eight and two locations, respectively.
To maximize model performance accuracy, a no-flow policy (water retention and release) was scheduled for the control gate at the Lotus-Outlet (SWL2) during the entire validation period. In fact, this particular flow policy was intended to improve the scheduling and communication of special events on the Lotus-Outlet (SWL2) and to maintain the primary goal of flood control and public safety.
Assessment of model performance
In this study, both qualitative (based on visual graphical technique) and quantitative (based on statistical measures) methods of assessment were combined but with greater emphasis placed on the statistical measures. The statistical measures for each point during the calibration and validation periods were evaluated (Tables 4 and 5, respectively). Furthermore, the predictive accuracy was tested using Pearson's distribution index (r2). The MIKE SHE simulations were considered optimal if the MAE, R, RMSE and Nash–Sutcliffe coefficient of efficiency (CE) were close to the values of 0, 1, 0 and 1, respectively. Obviously, the PIW model scored the highest values of R (0.765 − 0.927) and CE (0.748 − 0.828) during validation. However, large RMSE values were calculated for the flow rate during calibration at SWL2 (outlet; 0.766) and Langat River during validation (0.780). This bias was attributed to low or occasional absence of variation in the historical time-series datasets necessary for the simulation process. In contrast, the calibration period was characterized by some missing values in the observational data, which were therefore excluded from the statistical assessment accordingly. Actually, in terms of statistical evaluation and visual flow dynamic presentation, the best performance was achieved during the validation period, and for the surface water level and channel flow simulations in particular, rather than for the groundwater table. Further mathematical explanations of these statistics can be found in the literature (Gutpa et al. 1998; Legates & McCabe 1999).
RESULTS AND DISCUSSION
Owing to the unavailability of such information regarding actual past control strategies (e.g., the setting of gate opening and maximum rate of exchange) and exact operational schedules, the model could not represent such events accurately. This is illustrated in the calibrated water level hydrographs of the Lotus-Outlet (SWL2) (Figure 4(l)) and the closest lakes, i.e., Chalet, Typha and Lotus (Figure 4(f), 4(h) and 4(j)). These results reveal that drought conditions have caused the water level to drop by 20–60 cm.
The observed and simulated surface water level hydrographs during the validation period are illustrated in Figure 5. The results clearly reveal that the model response was much better during the validation period than during the calibration period in that the dynamics of the water level fluctuations are represented well. Furthermore, most of the simulated flow fluctuations match their observed counterparts well, especially for the upper and middle lakes, namely Visitor, Main, Tin, Crocodile and Hippo lakes (Figure 5(a)–5(e); Table 4). However, at the lower lakes (namely, Chalet, Typha and Lotus), it appears that the model missed capturing two anomalous peaks that occurred during two periods: (I) 4–11 January 2008 and (II) 22–25 March 2008 (Figure 5(f)–5(h)). These anomalies are justified by the occurrence of an event comprising four successive storms during the first period followed by another event comprising two storms during the second period, which accounted for 8.2 and 3.2%, respectively, of the total rainfall during the entire validation period. In fact, maintaining a no-flow policy at control gate (Lotus-Outlet) during the entire validation period contributed effectively to the satisfactory representation of the dynamic characteristics of the surface water flow in comparison with the calibration period. Thus, the overall slight underestimation, especially for the lower lakes, was mainly related to the accumulation of water in Lotus Lake after the storm events, which tended to exceed the discharge intake capacity of the broad-crested weir type of control gate.
Results of channel flow simulation during both periods of calibration and validation are presented in Tables 4 and 5, respectively, as well as Figure 6. There is a satisfactory agreement between the modeled flow and the observed hydrograph of the Inlet during the calibration period (Figure 6(a)), supported by a reasonably strong relationship based on the skewness index (r2) of 0.78, and by the very significant (P ≤ 0.01) prediction variables (Table 3).
These results, on the one hand, demonstrate the strong correlation between observed and simulated discharge, indicating the satisfactory performance of the PIW model. On the other hand, they emphasize the increased conformity of the modeled outcomes compared to the observed data, thereby confirming that the Lotus-Outlet displays a discrepancy due to the no-flow policy rather than the model simulation (Figure 6(b)), as evidenced by poor evaluation criteria (R = 0.45; CE = 0.42), values of distribution index (r2 = 021), as well as both of predictive variables of 0.35 and 0.27 for the slope and intercept, respectively (Table 4). This disagreement is attributed mainly to unscheduled operations that occurred at the weir of the Lotus-Outlet. Similar to the surface water level at this spot (Figure 4(l)), these operations obviously influenced the flow rate during both dry and wet periods, which accounted for the occurrence of many uncaptured limbs of the model-simulated hydrograph (Figure 6(b)). Because vandalism occurred for the automatic data logger at the Inlet spot, the reach of Langat River within the PIW catchment had been chosen in order to compensate the missing flow rate data during the validation period. The flow rate simulation result revealed a uniquely identical matching between measured and modeled values for the reach of Langat River stamped with some slight overestimation on some occasions (Figure 6(c)). Similarly, both patterns of modeled and measured flow rate at the Lotus-Outlet (Figure 6(d)) match identically each other, clear from remarkably (P < 0.01) high evaluation and prediction criteria of R = 0.927; CE = 0.828; r2 = 0.86 (Table 5). Noticeably, there is some logically slight misrepresentation at a few simulation points observed in the Lotus-Outlet hydrograph (Figure 6(d)) which was attributed to standardization of the dimensionless roughness value of ‘10’ to be assigned to the whole channels of the PIW. In fact, this assigned roughness value is considered slightly high relative to tropical climate conditions of the study area, where the plant-hold-back rivers are the common feature associated with the watercourses. This practice was adopted in order to overcome the unsteadiness condition of the PIW model during the calibration process that popped up because of the over-parameterization issue.
To further investigate the dynamic characteristics of flow at this problematic stimulation spot of Lotus-Outlet, the behavior of flow rate relative to rainfall events is illustrated in Figure 7. In response to rainfall events, the PIW model, therefore, succeeded to simulate the hydrological dynamics of the flow rate which included the overland flow as well specially from spots other than peat swamps where a very low roughness coefficient of ‘5’ was assigned in consistence with occurrence of many objects obstructing the flow dynamic including grassy vegetation and rough ground surface. Unlike small rainfall events, the model simulated the significant stormy events at its best satisfactory performance because of the insignificant contribution from OL; as there is no runoff over a thirst/unsaturated soil. Once again, the performance harmony of the validated PIW model has been interrupted by the over-accumulation of water in Lotus Lake after the storm events (Figures 6(d) and 7).
The groundwater head at BH3 was expected to have limited episodic influence in response to subsurface leakage from the nearby Main Lake, taking into account the substantial difference of groundwater head fluctuations at BH3 and those shown for this lake (Figures 4(e), 5(b) and 8(c)). Thus, no clear relation was evident between the groundwater head around BH3 and the water level of Main Lake that could be considered as aquifer–lake interaction in this part of the catchment. However, the downstream part of the catchment that extends from Lotus Lake to the reaches of the Langat River lies within the zone of influence of the Megasteel pumping well field, which strongly controls the groundwater heads.
The validated model showed satisfactory spatially distributed predictions of the dynamics of groundwater head with a performance somewhat similar to that achieved during calibration for BH3 (Figure 8(j); Tables 4 and 5). However, the nearly flat representation of the groundwater level in the simulation hydrograph at BH5 (Figure 8(i)) might reflect the occurrence of some expected uncertainties in this part of the modeled catchment, mostly associated with extensive groundwater withdrawal at the Megasteel well field. This assumption is strongly supported by the fact that the simulated groundwater table tended to rise by approximately 0.2 m at BH3 in comparison with the calibration period for the same piezometer, whereas the groundwater head at BH5 dropped by approximately 0.75 m during both wet (November–January) and dry (May–July) periods (Figure 8(e) and 8(i)).
CONCLUSIONS
The approach of multi-criteria simulation was applied to overcome limitations in both continuity and availability of the observational datasets used to model a complex data-limited wetland catchment. In fact, the hydrological connectivity of PIW (i.e., site specification) facilitates bridging gaps in historical datasets as the model calculates simultaneously all hydrological processes of different calibration/validation targets. Three calibration targets including surface water level, channel flow rate and groundwater head were specified for the period from 1 July 1999 to 31 October 2004. To maintain simulation of at least a single complete hydrological cycle, the calibrated coupled model was validated for the period from 1 August 2007–1 August 2008 using an accurate field dataset. The coupled model demonstrated satisfactory simulation performance for all the hydrodynamic components of the water cycle throughout the PIW catchment. The accuracy of the model applicability to PIW was assessed by high scores of the R and CE measures. Nonetheless, large RMSE values were calculated for channel flow rate during calibration at SWL2 (outlet) and during validation at Langat River (replacement of SWL1 inlet). In fact, this bias was an inevitable outcome following vandalism and occasional absence of variation in the historical time-series datasets necessary for the timely simulation process at SWL1 and SWL2, respectively.
Visual assessment and multi-criteria evaluation showed that both the calibrated and the validated models performed satisfactorily. However, the hydrographic visual assessment revealed that the dynamic characteristics, especially for surface water, were represented better by the validated model than by the calibrated model. The latter was noticeably influenced by unscheduled flow operation of the control gate at the Lotus-Outlet. Therefore, in this study, validation is defined only as substantiation that a site-specific model can perform simulations at a satisfactory level of accuracy; hence, no universal validity of the general model code is tested nor claimed. During validation, the groundwater dynamics in the area between Lotus Lake and the Langat River were influenced remarkably by groundwater abstraction at the Megasteel well field. The simulated channel flow for both the calibration and the validation period responded reasonably well to the input rainfall data. This showed that the water level in the lakes and channel system fluctuates in response to climatic variability. Moreover, because of the flat topography, it also has the ability to retain a significant proportion of the runoff, especially during storm events. Thus, the rate of surface water flow varies greatly from year to year and from event to event.
The flexibility of the MIKE SHE operating structure allows the use of any number of components depending on the availability of input data. Discontinuity of observational time-series datasets, the absence of records of actual historical and present flow control strategies, as well as the operational schedule, appear to control model accuracy. However, this modeling practice demonstrates that calibrating the model against more than two targets (e.g., surface water level, streamflow, runoff, groundwater heads and groundwater discharge) could improve the predictive accuracy of model performance. Unlike conventional simplified watershed models that cannot represent all the components of the water cycle, the unique feature of the MIKE SHE hydrological components is the integration of various hydrological processes at different timescales. Thus, the nature of the integration and the ability to account for both surface and subsurface flow systems and their interactions make the MIKE SHE model well suited for wetland studies.
ACKNOWLEDGEMENTS
The authors thank the University of Malaya (UM) and the Department of Irrigation and Drainage of Malaysia (DID).
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
CONFLICT OF INTEREST
The authors declare there is no conflict.