1. Introduction

In recent decades, climate change and socio-economic growth have had impacts on available water resources and food security in most parts of the world (Biemans et al. 2013; Lee, Bae 2015). Thus, more sophisticated management of existing water resources (Mereu et al. 2016; Wałęga et al. 2020), environments (Alam et al. 2018), and ecosystems (Yaegashi et al. 2018) is required to adapt to the impacts. A water reservoir’s operation should follow a rational policy to ensure adequate water provision for its various intended purposes without adverse effects. Figure 1 shows a conceptual diagram of a typical reservoir operation system in the context of feedback control theory. Hydrologic inputs, including precipitation, evapotranspiration, and runoff from the catchment, affect the reservoir’s water balance and the water supply to the command area (i.e. the downstream area serviced by the reservoir’s functions). The release from the reservoir as a control variable also contributes to the reservoir’s water balance and the water supply to the command area. The operator’s role is to implement a policy determining the release discharges based on the reservoir’s water level and the water demand from the command area. However, a more sophisticated system of reservoir operation might also account for hydrologic inputs and/or possible conveyance losses as parts of the information on which the policy is based.

Fig. 1.

Conceptual block diagram of a reservoir operation system.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g001_min.jpg

Labadie (2004), Rani and Moreira (2010), and Ahmad et al. (2014) reviewed a variety of methods for computationally determining suitable reservoir operation policies. Stochastic dynamic programming (SDP) theory is the general framework to establish reservoir operation policies, as explored decades ago (Heidari et al. 1971; Tejada-Guibert et al. 1993, 1995; Yakowitz 1982). The SDP typically provides optimal operation policies with thresholds of opening valves or switching pumps (Nop et al. 2021; Unami, Mohawesh 2018; Unami et al. 2019b). However, reservoir operators generally do not accept such theoretically generated policies but instead rely on their empirical knowledge to make decisions when unexpected risks are involved (El-Shafie et al. 2014). Therefore, attention should be paid to the actual policies that may be identified from historical data (Giuliani et al. 2014). Turner et al. (2021) firstly attempted to inventory reservoir operation policies in the United States.

This study aims at establishing a new approach to identifying nominal release policies implemented in a multi-purpose water reservoir. Bukit Merah Reservoir (BMR), located in the northern part of Perak State, Malaysia, was chosen as a study site. The water levels and release discharges of BMR were recorded daily from 2000 through 2011. Generalized additive models (GAMs; Chen, Tsay 1993) are applied for representing the operational policies for release discharges, which are assumed to have an annual period based on the day of the year and the water levels (Unami et al. 2019a). The backfitting algorithm works well to identify each contributing function of the GAMs (Breiman, Friedman 1985). However, the release policies identified with that method include spurious oscillations, which are removed with total variation regularization (TVR), known as the Rudin–Osher–Fatemi model (ROF) for the initials of its three authors (Rudin et al. 1992). The ROF model contains a scale parameter dominating the discrepancy of the reconstructed data from the original data that includes spurious oscillations (Osher et al. 2005). Further application of TVR in agricultural water management was discussed in Fadhil and Unami (2021). The identified release policies, which are considered nominal because they have been smoothed and regularized, are utilized to examine shifts in the operation of BMR during the period. Several factors stemming from climate change and socio-economic growth are inferred to have burdened the operation of BMR with more demanding release policies.

2. Materials and methods

2.1. Study site and datasets

Figure 2 shows the topography of BMR’s catchment and command areas, extracted from the Shuttle Radar Topography Mission (SRTM) digital elevation data (Farr et al. 2007), with the boundaries of sub-catchments demarcated with white lines. The climatic zone is tropical rainforest with a mean annual rainfall of more than 3000 mm. The bimodal annual rainfall pattern enables paddy rice cultivation twice a year. BMR has a catchment area of 480 km2; it supplies irrigation water to the paddy fields of the Kerian Irrigation Scheme (KIS), covering an area of 236 km2, and urban water for municipal and industrial purposes to the Kerian, Larut, Matang, and Selama Districts. The catchment includes four river systems: Merah, Jelutong, Selarong, and Kurau Rivers (Ismail, Najib 2011), located between 04 51 N and 05 10 N latitude and 100 38 E to 101 00 E longitude. The region is primarily rural with numerous riverine villages established along the middle and lower reaches of the rivers. The land use in the areas with elevations lower than 50 m is predominantly for tree plantations such as oil palm, rubber, and coconut. At the same time, the surrounding steep mountain slopes are covered with rainforests. The rainfall characteristics in the catchment area have been modeled with a stochastic generator (Fadhil et al. 2017). Analysis using rainfall-runoff models with fractional derivatives has shown that the hydrologic response in the catchment area of the Kurau River is unstable (Unami et al. 2021). The KIS’s 236 km2 command area consists of low-lying lands between the Strait of Malacca and BMR. Two primary irrigation canals, the Main Canal and the Selinsing Canal, convey the irrigation water by gravity. The urban water is taken from the Selinsing Canal at a pumping station located 6.5 km downstream of BMR. Dor et al. (2011) conducted a detailed hydro-geological study on the Selinsing Canal. Approximately 61% of water consumed in the paddy fields of the KIS originates from BMR, and the rest is from rainfall (Hamidon et al. 2015). Urban demand depends on the local population and its economy. The Department of Irrigation and Drainage (DID, 2011) projected that the growth rates of the population and economy in the state of Perak would decrease but remain positive until the year 2050. Intensive research on the operation of BMR has been conducted in the context of future climate change and SDP (Fadhil 2018). Besides the irrigation and urban demands, the operation of BMR needs to account for flood control and environmental hazards. There are two gated spillways of Ogee weir type with the same crest level, draining water from BMR to the original water course of the Kurau River, which meanders through the KIS with widths of 30-50 m.

Fig. 2.

The topography of the catchment and command areas of BMR, located in Peninsular Malaysia.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g002_min.jpg

The key dimensional parameters of BMR are summarized in Table 1. The dam embankment is at an elevation above the sea level (EL) of 11.28 m. The datum level is taken as the lowest level of the reservoir. The water level above the sea level is denoted by h (m). The storage volume of the reservoir when the water level is equal to h is denoted by V(h). The maximum capacity of the spillways at MWL is 565 m3/s. However, the bed slope of the Kurau River receiving the water from the spillways is about 1/10,000, implying insufficient discharge capacity.

Table 1

Key dimensional parameters of BMR.

ParameterDescription
Maximum water level (MWL) H (m)EL 10.25
Flood control level (FCL) ηF (m)EL 9.75
Normal water level (NWL) ηN (m)EL 8.69
Spillway crest level (SCL) ηS (m)EL 8.14
Intake gate crest level (IGL) ηG (m)EL 6.10
Datum level η0 (m)EL 2.00
Storage volume (106 m3)V(h) = 0.190 (hη0)3.17

A daily log of the operation of BMR in terms of water levels, rainfall depths, and release discharges into the two primary irrigation canals for the period from 2000 through 2011 was provided from the DID. The complete source data are available in a supporting information file.

Hamidon et al. (2015) estimated the irrigation demand Qirr in the KIS for each month. The urban demand Qurb was assumed to be constant at 1.70 m3/s (Anwar 2010). Then, the sum of the irrigation demand Qirr and the constant urban demand Qurb becomes the total demand discharge QD for each month. Table 2 summarizes the average monthly rainfall depths R during 2000-2011, the values of Qirr, Qurb, and QD, and the observed average monthly release discharge QR. However, it is noteworthy that the actual irrigation demand varies on a daily scale, according to the GIS-based estimation by Rowshon et al. (2003a-b).

Table 2

Monthly parameter values relevant to the water balance of BMR.

ParameterJanFebMarAprMayJunJulAugSepOctNovDec
Rainfall depth R (mm)243.0174.8294.3305.5155.0163.1170.3215.6296.9406.3356.6266.8
Irrigation demand Qirr (m3/s)11.350.3913.9915.6314.9513.6310.120.3511.5412.2312.5412.84
Urban demand Qurb (m3/s)1.701.701.701.701.701.701.701.701.701.701.701.70
Total demand QD (m3/s)13.052.0915.6917.3316.6515.3311.822.0513.2413.9314.2414.54
Observed average release discharge QR (m3/s)5.744.1824.0728.9127.1813.864.062.7923.3229.8828.5214.92

2.2. Generalized additive models for release discharges

The water level of BMR at the beginning of day t is denoted by ht. The release discharge into the primary irrigation canal k on the day t of the year y is denoted by Qy,tk. Here, k = 0 and k = 1 indicate the Main Canal and the Selinsing Canal, respectively. It is assumed that the operator’s nominal policy uk (t, ht) for release discharge into canal k, based on the information of t and ht, is yearly invariant. We introduce an additive decomposition structure into uk (t, ht) to write it as a GAM:

(1)
uk(t,ht)=ftimek(t)+fWLk(ht)+Ck

where ftimek and fWLk are functions that may be nonlinear, and Ck is a constant. Representing the day t as:

(2)
t=2πThedayoftheyear1Thenumberofdaysintheyear,

the function ftimek is assumed to be 2π-periodic. The objective here is to determine the functions ftimek and fWLk minimizing the errors:

(3)
εy,tk=Qy,tkut(t,ht)

in the sense of the least squares as well as the TV.

2.3. Backfitting algorithm

Firstly, the backfitting algorithm is applied for finding the functions to minimize the expected square error E[εt,yk]2 for all the observed data. The algorithm iteratively updates smoothing estimates f^timek and f^WLk of ftimek and fWLk, respectively. Histograms ϕtimek and ϕWLk are defined as:

(4)
ϕtimek(l)=1|Ωt(l)|tΩt(l)f^timek(t),Ωt(l)={t|t(lΔtΔt2,lΔt+Δt2)}

and

(5)
ϕWLk(l)=1|Ωh(l)|htΩh(l)f^WLk(ht),Ωh(l)={ht|ht(lΔhΔh2,lΔh+Δh2)}

for integers l with the increments Δt and Δh, respectively. Then, with initial guesses f^timek=f^WLk=0 for all t and ht, the histograms are utilized to update the smoothing estimates as:

(6)
f^timek(t)=Qy,tkQWLk(l)Ck,t(lΔtΔt2,lΔt+Δt2]

and

(7)
f^WLk(ht)=Qy,tkϕtimek(l)Ck,ht(lΔhΔh2,lΔh+Δh2]

where Ck are taken as the averages of all Qy,tk. This algorithm is known to converge under appropriate conditions.

2.4. TVR

After the backfitting algorithm has converged, the TVR is applied to each f^timek and f^WLk=0 to remove spurious oscillations. The TVR operated for a generic oscillating function f over a domain Ω of x minimizes the functional J of u:

(8)
J=Ω|u|dΩ+λ2Ω|uf|2dΩ

where λ is the scale parameter of the ROF model, and ∇ = ∂/∂x. The Euler-Lagrange equation for minimization of J in (8) is formally written as:

(9)
(u|u|)+λ(fu)=0.

We attempt to obtain an approximate solution to (9) by numerically solving the initial value problem of the singular diffusion equation:

(10)
ut=(u|u|)+λ(fu),u=fatt=0

It is assumed that the values of f are determined at n different points of Ω. Let those points be xi for i = 0, 1, …, n – 1 such that x0 < x1 < … < xn–1. Applying the standard Galerkin finite element scheme with piecewise linear bases to (10) results in the initial value problem of the ordinary differential equation:

(11)
dudt=M1Aσ+λ(fu),u=fatt=0

where u ∈ ℝn is the vector whose ith entry is uι, which is the approximation of u at xi, f ∈ ℝn is the vector whose ith entry is f at xi, σ ∈ ℝn−1 is the vector whose ith entry is the sign of ui+1ui, M ∈ ℝn×n is the mass matrix, and A ∈ ℝn×(n−1) is the stiffness matrix. Two types of boundary conditions are considered: periodic and Neumann. For the periodic boundary condition, the matrices become:

(12)
M=(x1xn13x1x0600x0xn16x1x06x2x03x2x16000xixi16xi+1xi13xi+1xi6000xn2xn36xn1xn33xn1xn26x0xn1600xn1xn26x0xn23)

and

(13)
A=(1001100110011001)

For the Neumann boundary condition, the matrices become:

(14)
M=(x1x06x1x06000x1x06x2x03x2x16000xixi16xi+1xi13xi+1xi6000xn2xn36xn1xn33xn1xn26000xn1xn26xn1xn26)

and

(15)
A=(1000100110010001)

3. Results

3.1. Identification of nominal release policies and reconstruction of release discharges

Two five-year subperiods, 2001-2005 and 2007-2011, were extracted from the whole period 2000-2011. Consecutive use of the backfitting algorithm and the TVR identified nominal release policies from the observed release discharges during each of the two subperiods. The increments are chosen as Δt=2π365.25 (rad) and Δh = 0.01 (m). Five cases of the scale parameter λ = 10i/2 (i = 0, 1, 2, 3, 4) are considered. The numerical solution of the initial value problem (10) successfully converges to a steady state for each case of λ. Figures 3 and 4 compare the processes of TVR applied to ftime0 for the Main Canal during the two subperiods 2001-2005 and 2007-2011. Figures 5 and 6 compare the processes of TVR applied to fWL0 for the Main Canal during the two subperiods. The value of C0, which is equal to the average of all release discharges to the Main Canal, is 12.01 and 15.60 for the subperiods 2001-2005 and 2007-2011, respectively. Then, release discharges to the Main Canal are reconstructed with the nominal release policies identified from the observed data during the two subperiods, as shown in Figures 7 and 8. For example, with λ = 10–2, Figures 3 and 5 indicate that the values of ftime0 and fWL0 are –8.29 and –0.28, respectively, if the water level is at SCL (EL 8.14) on January 1st during the subperiod 2001-2005. In that case, the nominal release discharge becomes equal to ftime0 + fWL0 + C0 = −8.29 − 0.28 + 12.01 = 3.44 (m3/s). Figures 9-13 show the results similarly computed for the Selinsing Canal as in Figure 3 through Figure 8 for the Main Canal. The value of C1, which is equal to the average of all release discharges to the Selinsing Canal, is 6.32 and 6.86 for the subperiods 2001-2005 and 2007-2011, respectively.

Fig. 3.

The process of TVR applied to the function of time for the Main Canal during the subperiod 2001-2005.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g003_min.jpg
Fig. 4.

The process of TVR applied to the function of time for the Main Canal during the subperiod 2007-2011.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g004_min.jpg
Fig. 5.

The process of TVR applied to the function of WL for the Main Canal during the subperiod 2001-2005.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g005_min.jpg
Fig. 6.

The process of TVR applied to the function of WL for the Main Canal during the subperiod 2007-2011.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g006_min.jpg
Fig. 7.

Reconstructed release discharges to the Main Canal with the 2001-2005 policies.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g007_min.jpg
Fig. 8.

Reconstructed release discharges to the Main Canal with the 2007-2011 policies.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g008_min.jpg
Fig. 9.

The process of TVR applied to the function of time for the Selinsing Canal during the subperiod 2001-2005.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g009_min.jpg
Fig. 10.

The process of TVR applied to the function of time for the Selinsing Canal during the subperiod 2007-2011.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g010_min.jpg
Fig. 11.

The process of TVR applied to the function of WL for the Selinsing Canal during the subperiod 2001-2005.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g011_min.jpg
Fig. 12.

The process of TVR applied to the function of WL for the Selinsing Canal during the subperiod 2007-2011.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g012_min.jpg
Fig. 13.

Reconstructed release discharges to the Selinsing Canal with the 2001-2005 policies.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g013_min.jpg
Fig. 14.

Reconstructed release discharges to the Selinsing Canal with the 2007-2011 policies.

https://www.mhwm.pl/f/fulltexts/154592/MHWM-9-0010-g014_min.jpg

3.2. Choice of the scale parameter

The daily oscillations in f^timek and f^WLk are successfully removed with TVR, while the discrepancy of ftimek and fWLk from them becomes large as the scale parameter λ becomes small. The TV in the reconstructed release discharges tends to decrease for smaller λ. Under the assumption that the release policies are based on the information of t and ht, the nominal release policies identified by the large-scale parameter λ = 1, are well representative of t and ht, and thus, considered the most reasonable and realistic.

4. Discussions

4.1. Comparison of the two subperiods

The comparison between the observed and the reconstructed release discharges during the whole period validates the release policies identified from the data during each of the two subperiods. Significant increases in the peaks and the TV of the observed release discharges in both canals can be seen in 2009-2011. Consequently, the reconstructed release discharges with the 2001-2005 policies barely approximate the observed ones in 2009-2011. This indicates that the operator changed the release policies from the subperiod 2001-2005 to the subperiod 2007-2011.

Referring again to Figure 1, the variable hydrological input, which is under the considerable influence of climate change, is considered the primary cause of changes in the release policies controlling the water balance in the reservoir to satisfy the water demand.

The changes in the functions ftimek from the subperiod 2001-2005 to the subperiod 2007-2011 imply the shifts in the actual irrigation demand’s annual patterns, which are susceptible to rapidly developing agronomic practices such as the choice of varieties to adapt the climate change and to meet the increasing food demand. In the Main canal (k = 0), the two peaks of early-March through early-May and early-September through late-October shifted to mid-March through late-May and mid-September through early-December, respectively. In the Selinsing canal (k = 1), the two peaks of early May through late June and late October through late December shifted to mid-April through mid-June and early-November through mid-January, respectively. As a result, May became the most critical month in the subperiod 2007-2011, having the least rainfall and the overwrapping peaks of the release discharges in the two primary irrigation canals. This increasing criticality of May in the water balance of BMR supports the simulation-based prediction by Hamidon et al. (2015). Furthermore, it is noteworthy that the change in rainfall patterns between the two subperiods in the catchment area of BMR is just a symptom of catastrophic regime shifts in the coming decades, according to the recent climatological studies using GCMs (Adib et al. 2020; Adib, Harun 2022; Adib et al. 2022).

The function fWL0 for the Main canal is mostly monotone increasing with respect to ht, and its increment was amplified from the subperiod 2001-2005 to the subperiod 2007-2011. The function fWL1 for the Selinsing canal was small for the higher water levels during the subperiod 2001-2005 but became not so much during the subperiod 2007-2011. In other words, the two primary irrigation canals became more utilized for draining excess water from BMR in the subperiod 2007-2011 than in the subperiod 2001-2005. This allocation of the two primary irrigation canals for the drainage purpose might be attributed to the small capacity of Kurau River, which cannot manage the discharges required for reasonable flood control.

4.2. Comparison with other approaches

The first key feature of the approach developed here is to decompose a release policy into possibly singular functions of a single independent variable with the additive structure, unlike the standard procedure identifying release policies as quite regular functions of multiple independent variables (Turner et al. 2021). As discussed in the previous subsection, the decomposition of release policies illustrates the two aspects of the irrigation demand’s annual patterns and the hydraulic structures’ functions.

The second key feature is the TVR applied to the decomposed possibly singular functions, breaking through a widespread obsession that policy identification is “to minimize some distance metric between historical releases and modeled ones,” as in Giuliani et al. (2014). This study is the first to succeed in quantitatively treating the indecisive attitude of a reservoir operator as the spurious oscillations to be removed with the TVR. The release policies extracted with the TVRs are the nominal parts. Overgaard (2019) provides mathematical proof to endorse the ROF model in one dimension used in this study for TVR. However, the lack of an established method to choose the best scale parameter λ is a practical limitation of this approach. A follow-up study shall be conducted to validate the λ values with the empirical knowledge of the reservoir operator.

A sophisticated approach is being developed to combine the backfitting algorithm with TVR (Yang, Tan 2018). Nevertheless, our approach, opting for the sequential use of the two procedures, is advantageous because it is not computationally demanding.

5. Conclusions

The consecutive use of the backfitting algorithm and TVR successfully identified the nominal release policies implemented in BMR during each of the two subperiods as the GAMs. The additive structure decomposed the varying parts of the nominal release policies into the two functions. The scale parameter of the TVR controlling the removal of spurious oscillations should be chosen so that the nominal release policies are sensitive to the time and water levels. The results have been utilized to compare the two subperiods, indicating the significant changes in the operator’s release policies. The rapidly developing agronomic practices and the insufficient discharge capacity of the Kurau River receiving the water released from the spillways might be seriously burdening the operation of BMR. We have also addressed the situation in another study (Fadhil 2018), so that the reservoir operator is aware of future hydro-meteorological information, and deducing optimal release policies of BMR based on the SDP theory.