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.

  • 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.

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

The region of the Paya Indah Wetlands (PIW) is located in the Kuala Langat District in the State of Selangor, Malaysia (Figure 1). It covers an area of ∼242 km2 and it encompasses myriad ecosystems, e.g., old tin mining ponds and peat swamp forest. The study area forms part of the Langat River Basin, which consists of metamorphosed sandstone, shale, mudstone and schist (Omar et al. 1999; Minerals and Geoscience Department of Malaysia 2002). In the low flatlands, thick Quaternary deposits that lie on the bedrock are distributed along the coastal fringe.
Figure 1

Detailed location map of Paya Indah Wetland (PIW) Catchment.

Figure 1

Detailed location map of Paya Indah Wetland (PIW) Catchment.

Close modal

Input datasets

Digitized topography dataset (Jurukur Permata Malaysia 2003) was used to generate a 20-m digital elevation model (DEM) for PIW (Figure 2(a)). Data for soil map were based on Wong (1966) (Figure 2(b); CorelDRAW 2020; ESRI 2020), whereas the spatial hydraulic properties datasets including hydraulic conductivity (Ks) for each soil type (namely, Peat (6.43 × 10−5 m/s), Mined and soil association (4.42 × 10−8 m/s), Serdang-Bungor Munchong series (4.53 × 10−7 m/s), Selangor-Kanchung series (3.95 × 10−7 m/s) and Perang series (6.63E × 10−7 m/s)) were measured in the field. Rainfall time-series sets were from October 1999 to September 2004 (calibration) and August 2007 to August 2008 (validation). The rainfall input for the model was distributed spatially according to a weighted method in which the total rainfall was calculated from the measured rainfall and area-weighted factors (Figure 2(c); Table 1). The data were then converted into the MIKE SHE time-series format.
Table 1

Rainfall area-weighted factors for the Paya Indah Wetland Catchment (area 242.21 km2)

StationWeighted factor (%100)
Sungai Manggis 40.10 
Bukit Cheeding 36.17 
Perang Besar 20.78 
KLIA 2.95 
StationWeighted factor (%100)
Sungai Manggis 40.10 
Bukit Cheeding 36.17 
Perang Besar 20.78 
KLIA 2.95 
Figure 2

(a) Distribution of spot level dataset and topographic model of PIW catchment. (b) Soil sampling and in-situ measurements spots and modeled soil map for PIW catchment. (c) Estimated rainfall fields for the catchment using the Thiessen polygon method. (d) Land use map of the Paya Indah Wetland Catchment. (e) Flooded areas within the Paya Indah Wetland Catchment: Lakes names and their representing code number.

Figure 2

(a) Distribution of spot level dataset and topographic model of PIW catchment. (b) Soil sampling and in-situ measurements spots and modeled soil map for PIW catchment. (c) Estimated rainfall fields for the catchment using the Thiessen polygon method. (d) Land use map of the Paya Indah Wetland Catchment. (e) Flooded areas within the Paya Indah Wetland Catchment: Lakes names and their representing code number.

Close modal

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.

Table 2

Properties of vegetation within the Paya Indah Wetland Catchment

Vegetation developmentTownship
Citrus
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 200 0.50 4.5 1,250 
Harvest 365 200 0.50 365 4.5 1,250 
Vegetation developmentSugarcane
Pasture
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 0.0 1.0 500 600 
Early season 60 2.0 1,000 150 600 
Mid season 90 2.5 1,500 240 600 
Late season 120 3.5 1,500 330 600 
1st harvest 150 4.5 1,500 365 600 
2nd harvest 210 5.5 1,500     
3rd harvest 365 6.0 1,500     
Vegetation developmentTruck crop
Grass
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 4.5 750 750 
Harvest 365 4.5 750 365 750 
Vegetation developmentShrub
Marsh
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 750 750 
Harvest 365 750 150 750 
     365 750 
Vegetation developmentSparse forest
Oil palm
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 750 2.5 0.85 
Harvest 365 750 365 2.5 0.90 
Vegetation developmentTownship
Citrus
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 200 0.50 4.5 1,250 
Harvest 365 200 0.50 365 4.5 1,250 
Vegetation developmentSugarcane
Pasture
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 0.0 1.0 500 600 
Early season 60 2.0 1,000 150 600 
Mid season 90 2.5 1,500 240 600 
Late season 120 3.5 1,500 330 600 
1st harvest 150 4.5 1,500 365 600 
2nd harvest 210 5.5 1,500     
3rd harvest 365 6.0 1,500     
Vegetation developmentTruck crop
Grass
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 4.5 750 750 
Harvest 365 4.5 750 365 750 
Vegetation developmentShrub
Marsh
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 750 750 
Harvest 365 750 150 750 
     365 750 
Vegetation developmentSparse forest
Oil palm
D.O.GLAIRDKcD.O.GLAIRDKc
Planting 750 2.5 0.85 
Harvest 365 750 365 2.5 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

A 2D conceptual geological model consisting of three layers was developed (Figure 3) based on the input data (Figure 2). The hydraulic properties of aquifers and the vadose zone and their spatial distribution beneath the PIW were determined from the field permeability test conducted within PIW catchment, and geological profiles based on previous comprehensive hydrogeological and soil investigations (Rahim et al. 2008, 2010) conducted in the adjacent Ampar Tenang area (Figure 1). In fact, over the previous few decades of extensive tin mining activities within the study area, the original geological formation, especially around the lake areas, has been covered completely by secondary mining deposits.
Figure 3

Conceptual model for Paya Indah Wetland (PIW).

Figure 3

Conceptual model for Paya Indah Wetland (PIW).

Close modal

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.

Water movement in the unsaturated zone

As a basic assumption, water flow in the unsaturated zone can occur as a Darcy flow within the soil matrix, as a gravity flow in distinct macropores (macropore flow), or as a matrix flow regime that can be described by Richards's equation:
(1)
where C is the soil water capacity (mm−1), Ψ is the pressure head (mm), K is the saturated hydraulic conductivity (mm s−1), Z is the gravitational head (mm) and S is the root extraction sink term (s−1).
Therefore, water movement in unsaturated soil is governed by an equation that can be derived by combining Darcy's law with the principle of mass conservation (Richards 1933). The pressure head form of the equation for a 1D vertical flow is:
(2)
where C(h) is the differential water capacity (∂θ/∂h) (L−1), θ is the volumetric water content (L3/L3), h is the soil water pressure head (matrix) (L), t is the time (T), z is the vertical coordinate (L), K is the isotropic hydraulic conductivity (L/T) and S is the sink term that represents root water extraction (L3/L3/T). In this study, the simplest approach was applied, which involved establishing a mass balance for the root zone and calculating the average moisture content for the entire root zone. If a sufficient amount of water is available within the root zone or in the capillary fringe, the simulated actual evapotranspiration is equal to the potential rate. In a water shortage situation, the potential rate reduces to a function of the available soil moisture within the root zone. Although the model has the benefits of being simple to use and fast to run, it is still based on physical soil properties. At this juncture, during the calibration process, manual adjustments were necessary to select the optimal values of some sensitive parameters, including soil material of vadose zone and aquifer properties. The lower boundary condition for the unsaturated zone in the model simulates the groundwater potential head. Hence, the depth of the unsaturated zone changes dynamically depending on climatic and hydrologic conditions.

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.

Table 3

Optimal values and ranges of sensitive parameters used for PIW model calibration

ModuleSensitive parameterOptimal values and ranges
Evapotranspiration Interception coefficient Cint 0.5 mm 
C0.3 
C0.2 
C30 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 
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) 
ModuleSensitive parameterOptimal values and ranges
Evapotranspiration Interception coefficient Cint 0.5 mm 
C0.3 
C0.2 
C30 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 
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) 
Table 4

Evaluation criteria for the calibrated model

Elements of multi-objective calibration strategy
Coordinates (Cassini system)
Evaluation criteriaa
PIW prediction variablesb
Calibration targetCategoryName of calibration pointXYMAERMSEcSTDresRCESlopeInterceptr2
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 targetCategoryName of calibration pointXYMAERMSEcSTDresRCESlopeInterceptr2
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).

Table 5

Evaluation criteria for the validated model

Elements of multi-objective validation strategy
Coordinates (Cassini system)
Evaluation criteriaa
PIW prediction variablesb
Validation targetCategoryName of validation pointXYMAERMSEcSTDresRCESlopeInterceptr2
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 targetCategoryName of validation pointXYMAERMSEcSTDresRCESlopeInterceptr2
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).

The model was calibrated and validated against the water level at 12 (Figure 4) and eight locations, respectively (Figure 5). Concurrently, the flow rate was included as an additional calibration target in the simulation process at three further locations to strengthen the predictive accuracy of the model (Figures 6 and 7). Satisfactory performance (Table 4) was obtained in nine calibration locations at which agreement between observed and simulated water levels was well matched (correlation (R) ≥ 0.75) during the entire calibration period. The observable overestimation and underestimation trends, which occurred only during short time intervals, were attributed to the normal uncertainty associated with modeling that could possibly be due to scaling problems (Legates & McCabe 1999). However, some sharp overestimated fluctuations in water level were simulated, causing considerable mismatch (0.45 ≥ R ≥ 0.65) between observed and simulated values during the following periods: (1) November − October 1999 (Figure 4(f), 4(h) and 4(j)), (2) November − March 2003 and (3) August − September 2004 (Figure 4(l)). It appears that an appropriate water level for the PIW lake system during special events (e.g., storms and droughts) could be maintained only if the control gate at the Lotus-Outlet (SWL2) was opened and closed as required by the Water Release and Reduced Flow Policy.
Figure 4

Hydrographs for simulation of surface water level during the calibration period.

Figure 4

Hydrographs for simulation of surface water level during the calibration period.

Close modal
Figure 5

Hydrographs for simulation of surface water level during the validation period.

Figure 5

Hydrographs for simulation of surface water level during the validation period.

Close modal
Figure 6

Simulation of PIW channel flow at the inlet and outlet spots during calibration and validation periods, respectively. Hydrographs (a) and (b) represent calibrated flow at SWL1 and SWL2, while (c) and (d) represent validated flow at Langat River and SWL2, respectively.

Figure 6

Simulation of PIW channel flow at the inlet and outlet spots during calibration and validation periods, respectively. Hydrographs (a) and (b) represent calibrated flow at SWL1 and SWL2, while (c) and (d) represent validated flow at Langat River and SWL2, respectively.

Close modal
Figure 7

Hydrodynamic fluctuations of the modeled flow rate at Lotus-Outlet relative to rainfall events.

Figure 7

Hydrodynamic fluctuations of the modeled flow rate at Lotus-Outlet relative to rainfall events.

Close modal

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 results for groundwater heads during the calibration and validation periods are shown in Figure 8 and Tables 4 and 5. Five out of eight locations show good agreement between simulated and observed groundwater heads during the calibration period in terms of dynamic fluctuation. It was found that the rise-and-fall patterns in the piezometers at BH1, BH2 and BH3 (Figure 8(a)–8(c), respectively) are much more dynamic in the upland areas in the northern upper part of the modeled catchment in comparison with BH6, BH7 and BH8 (Figure 8(f)–8(h), respectively) that show gradual fluctuation. It seems that the groundwater head in the upper part of the catchment rose and declined abruptly in response to the periodic flow exchange between the saturated and unsaturated zones caused by seepage from the overland flow and inflows from the surrounding Kuala Langat swamp forest. These fluctuations caused the groundwater table across the modeled catchment to drop between 0.45 and 0.65 m during the dry season. In fact, while the measurements are point values collected at groundwater piezometers, the model simulations are representative of average groundwater head within an area of ∼242 km2.
Figure 8

Simulation of groundwater head at PIW. Hydrographs (a)–(h) represent the calibration period. Hydrographs (i) and (j) represent the validation period.

Figure 8

Simulation of groundwater head at PIW. Hydrographs (a)–(h) represent the calibration period. Hydrographs (i) and (j) represent the validation period.

Close modal

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)).

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.

The authors thank the University of Malaya (UM) and the Department of Irrigation and Drainage of Malaysia (DID).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Adeyeri
O. E.
,
Laux
P.
,
Arnault
J.
,
Lawin
A. E.
&
Kunstmann
H.
2020
Conceptual hydrological model calibration using multi-objective optimization techniques over the transboundary Komadugu-Yobe basin, Lake Chad Area, West Africa
.
J. Hydrol.: Reg. Stud.
27
,
100655
.
Ahmadi
M.
,
Arabi
M.
,
AscoughII
J. C.
,
Fontan
D. G.
&
Engel
B. A.
2014
Toward improved calibration of watershed models: Multisite multiobjective measures of information
.
Environ. Modell. Software
59
,
135
145
.
Alminagorta
O.
,
Rosenberg
D. E.
&
Kettenring
K. M.
2016
Systems modeling to improve the hydroecological performance of diked wetlands
.
Water Resour. Res.
52
,
7070
7085
.
Beaulieu
J. J.
,
Balz
D. A.
,
Birchfield
M. K.
,
Harrison
J. A.
,
Nietch
C. T.
,
Platz
M. C.
,
Squier
W. C.
,
Waldo
S.
,
Walker
J. T.
,
White
K. M.
&
Young
J. L.
2018
Effects of an experimental water level drawdown on methane emissions from a eutrophic reservoir
.
Ecosystems
21
,
657
674
.
Boussinesq
J. M.
1904
Recherches théoriques sur l’écoulement des nappes d'eau infiltrées dans le sol et sur le modell des sources
.
J. Math. Pures Appl.
5
(
10
),
5
78
.
DHI
2004
MIKE SHE User Manual. Hørsholm, Denmark: Danish Hydraulic Institute. In Dingman, S.L. 2002 Physical Hydrology, 2nd edn. Upper Saddle River, NJ: Prentice Hall
.
ESRI
2020
ArcGIS Pro. Desktop 10.8. Environmental Systems Research Institute, Inc. Redlands. Available from: https://www.esri.com/en-us/arcgis/products/arcgis-desktop/resources.
FAO
2012
Crop evapotranspiration – Guidelines for computing crop water requirements – FAO Irrigation and Drainage Paper 56. Food and Agriculture Organization (Rome, 1998). Available from: http://www.fao.org/docrep/X0490E/X0490E00.htm.
Gutpa
H. V.
,
Sorooshian
S.
&
Yapo
P. O.
1998
Toward improved calibration of hydrologic models: multiple and noncommensurable measures of information
.
Water Resour. Res.
34
(
4
),
751
763
.
Jurukur Permata Malaysia
2003
Spot Levels Survey for Payah Indah. Reports # JPM (SEL) 630/2003/L1 to L10 and JPM(SEL) 360/2003/N1 to N13. Jurukur Permata Malaysia (Kuala Lumpur)
.
Liu
J.
,
Liu
T.
,
Bao
A.
,
De Maeyer
P.
,
Feng
X.
,
Miller
S. N.
&
Chen
X.
2016
Assessment of different modelling studies on the spatial hydrological processes in an arid alpine catchment
.
Water Resour. Manage.
30
,
1757
1770
.
Minerals and Geoscience Department of Malaysia
2002
The study on sustainable groundwater resources and environmental management for Langat Basin in Malaysia. Interim Report No. 1. (Kuala Lumpur)
.
Pechlivanidis
I. G.
,
Jackson
B. M.
,
Mcintyre
N. R.
&
Wheater
H. S.
2011
Catchment scale hydrological modelling: A review of model types, calibration approaches and uncertainty analysis methods in the context of recent developments in technology and applications
.
Global NEST J.
13
(
3
),
193
214
.
Prucha
B.
,
Graham
D.
,
Watson
M.
,
Avenant
M.
,
Esterhuyse
S.
,
Joubert
A.
,
Kemp
M.
,
King
J.
,
le Roux
P.
,
Redelinghuys
N.
,
Rossouw
L.
,
Rowntree
K.
,
Seaman
M.
,
Sokolic
F.
,
van Rensburg
L.
,
van der Waal
B.
,
van Tol
J.
&
Vos
T.
2016
MIKE-SHE integrated groundwater and surface water model used to simulate scenario hydrology for input to DRIFT-ARID: The Mokolo River case study
.
Water SA
42
(
3
),
384
398
.
Rahim
B. E. A.
,
Yusoff
I. B.
,
Samsudin
A. R.
,
Yaacob
W. Z.
&
Rafek
A. G.
2008
Heavy metal contamination of soil beneath a waste disposal site at Dengkil, Selangor, Malaysia
.
Soil Sediment Contam.
17
(
5
),
449
466
.
Rahim
B. E. A.
,
Yusoff
I. B.
,
Samsudin
A. R.
,
Yaacob
W. Z.
&
Rafek
A. G.
2010
Deterioration of groundwater quality in the vicinity of an active open-tipping site in West Malaysia
.
Hydrogeol. J.
18
(
4
),
997
1006
.
Sandu
M.-A.
&
Virsta
V.
2015
Applicability of MIKE SHE to simulate hydrology in Argesel River Catchment
.
Agric. Agric. Sci. Procedia
6
,
517
524
.
Vansteenkiste
T.
,
Tavakoli
M.
,
Ntegeka
V.
,
Willems
P.
,
De Smedt
F.
&
Batelaan
O.
2013
Climate change impact on river flows and catchment hydrology: A comparison of two spatially distributed models
.
Hydrol. Processes
27
,
3649
3662
.
Wong
I. F. T.
1966
Reconnaissance Soil Survey of Selangor
.
Ministry of Agriculture & Co-operatives
,
Kuala Lumpur
,
Malaysia
, pp.
11
28
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).