Development of the Indus River System Model to Evaluate Reservoir Sedimentation Impacts on Water Security in Pakistan

: Pakistan’s society and economy are highly dependent on the surface and groundwater resources of the Indus River basin. This paper describes the development and implementation of a daily Indus River System Model (IRSM) for the Pakistan Indus Basin Irrigation System (IBIS) to examine the potential impact of reservoir sedimentation on provincial water security. The model considers both the physical and management characteristics of the system. The model’s performance in replicating provincial allocation ratios is within 0.1% on average and the modeling of water ﬂow at barrages and delivered to irrigation canal commands is in agreement with recorded data (major barrage NSE 0.7). The average maximum volumetric error for the Tarbela and Mangla reservoirs are respectively 5.2% and 8.8% of mean annual inﬂow. The model showed that a 2.3 km 3 reduction in storage volume since 1990 equates to approximately 1.3 km 3 i.e., a 4–5% reduction in irrigation deliveries, respectively, for Punjab and Sindh in the dry (Rabi) season. This decline indicates that without further augmentation of system storage, the Rabi season supplies will continue to be further impacted in the future. This paper demonstrates the suitability of IRSM for exploring long term planning and operational rules and the associated impacts on water, food and energy security in Pakistan.


Introduction
Pakistan is facing severe challenges related to water, food and economic security. It relies utterly on the water of the River Indus and the associated groundwater systems. It operates the world's largest continuous irrigation system, the Indus Basin Irrigation System (IBIS), which supports food production, energy generation, as well as stock, domestic and industrial supply for the nation [1]. Predicting future seasonal inflows and the subsequent seasonal allocation, delivery and distribution of surface water resources throughout the IBIS is a laborious and complex process for water managers [2,3]. Water security issues in the basin include changing resource availability [4][5][6][7][8][9][10][11][12][13], growth in demand from competing users [14][15][16][17][18] and the loss of storage capacity through sedimentation [4,19]. The 2012 Friends of Democratic Pakistan Water Sector Task Force report [20] highlighted the need to build the knowledge and capacity to manage the current and future water security needs of Pakistan, with Young, et al. [21] highlighting the need to model the system accurately and transparently as a means to better manage water security risks. This study focusses on modelling both the physical and water-sharing arrangements of the IBIS to investigate the potential impacts on seasonal water supply security resulting from the loss of storage capacity through reservoir sedimentation.
Pakistan to ensure that irrigation areas in southern Punjab traditionally supplied by the eastern rivers could now be supplied from western rivers (Figure 1). the 1960 treaty precipitated construction of the two major reservoirs and numerous link canals in Pakistan to ensure that irrigation areas in southern Punjab traditionally supplied by the eastern rivers could now be supplied from western rivers (Figure 1).
For Pakistan, the Indus River and the associated groundwater system is the main source of water supply. Approximately 71% of the 167 km 3 (1976-2020) mean annual inflow to the IBIS is used for agriculture and to a lesser extent for industrial and domestic consumption. Surface water contributes approximately 60% of the irrigation water use in Pakistan [12], with the remainder being made up from groundwater, predominantly in Punjab and in isolated pockets of Sindh province [34], where groundwater supplies are fresh. Figure 1. The Indus River system, showing major reservoirs, headworks, provincial boundaries, the Indus Basin Irrigation System (IBIS) and the Indus and Jhelum-Chenab (J-C) management zones (systems) in Pakistan. Tarbela, Mangla, Marala, Nowshera, Balloki and Sulemanki are considered remote inflow measurement (rim) stations for IBIS. For a detailed schematic diagram of the IBIS (depicting the current storage and canal system as of 2020), refer to Figure A1 in Ahmad et al. [12].

Climate and Seasons
The climate of this region is derived from the South Asian monsoon and western disturbances (WDs) [35]. The monsoon occurs during the June to September period and is responsible for most of Pakistan's precipitation. This period also coincides with the period of maximum snow and glacier melting in the Upper Indus Basin (UIB).
Most of the agriculture is practiced in the lower Indus plains. Summer spans April through to September, with maximum temperature ranging from 21 °C to 49 °C. Winter extends from December to February, with maximum air temperature ranging from 25 °C to 27 °C during day, with colder overnight temperatures. The mean annual rainfall in this area is about 250 mm/year, of which about 75% occurs between June to September. There are two distinct cropping seasons in Pakistan-Kharif (April to September) and Rabi (October to March). The dominant crop in the dry (winter) Rabi season is wheat, whereas cotton, rice and maize are the predominant crops during the wet Kharif (summer) season. Sugarcane and different fodder crops are grown throughout the year [36].  Figure A1 in Ahmad et al. [12].
For Pakistan, the Indus River and the associated groundwater system is the main source of water supply. Approximately 71% of the 167 km 3 (1976-2020) mean annual inflow to the IBIS is used for agriculture and to a lesser extent for industrial and domestic consumption. Surface water contributes approximately 60% of the irrigation water use in Pakistan [12], with the remainder being made up from groundwater, predominantly in Punjab and in isolated pockets of Sindh province [34], where groundwater supplies are fresh.

Climate and Seasons
The climate of this region is derived from the South Asian monsoon and western disturbances (WDs) [35]. The monsoon occurs during the June to September period and is responsible for most of Pakistan's precipitation. This period also coincides with the period of maximum snow and glacier melting in the Upper Indus Basin (UIB).
Most of the agriculture is practiced in the lower Indus plains. Summer spans April through to September, with maximum temperature ranging from 21 • C to 49 • C. Winter extends from December to February, with maximum air temperature ranging from 25 • C to 27 • C during day, with colder overnight temperatures. The mean annual rainfall in this area is about 250 mm/year, of which about 75% occurs between June to September. There are two distinct cropping seasons in Pakistan-Kharif (April to September) and Rabi (October to March). The dominant crop in the dry (winter) Rabi season is wheat, whereas cotton, rice and maize are the predominant crops during the wet Kharif (summer) season. Sugarcane and different fodder crops are grown throughout the year [36].

IBIS Infrastructure
There are 6 main remote inflow measurement (rim) stations in the IBIS-Kabul River at Nowshera, Indus River at Tarbela, Jhelum River at Mangla, Chenab River at Marala, Ravi River at Balloki and Sutlej River at Sulemanki. The system's operation relies on two large supply reservoirs, Tarbela and Mangla, and a network of regulating infrastructure, consisting of Chashma reservoir with 0.34 km 3 storage and 15 river barrages ( Figure 1).
Tarbela reservoir is located on the Indus River and can store about 12% of the mean annual flow. This reservoir has been operational since 1976. Originally, it had 12 km 3 of live storage capacity, however after 44 years of operations and sedimentation, the live storage in 2020 has been reduced to 7.38 km 3 .
Mangla reservoir can store up to 36% of the mean annual flow of the Jhelum River. This reservoir has been operational since 1965 and the initial live storage of this reservoir was 7.25 km 3 in 2005, although due to sedimentation, the reservoir capacity was reduced to 5.77 km 3 . The reservoir was subsequently raised in 2009 to provide 9.12 km 3 storage. According to a recent (2018) available survey [37], the Mangla reservoir capacity is now at 9.08 km 3 .
River flows and reservoir releases are regulated to irrigation command areas through command canals that offtake from river barrages. Transfers between rivers also take place at barrages through link canals. From a water management perspective, the management system is broken into two sub-systems-Indus and Jhelum-Chenab (J-C).

Water Inflow Forecast, Allocation and Sharing Methodology
Due to variations in flow regime in the different tributaries of the Indus and limited storage capacity in Tarbela and Mangla reservoirs, a complex planning process is undertaken on a seasonal basis to equitably share the water. Figure 2 shows the mean 10-day inflow and demand patterns for (a) the Indus system and (b) Jhelum-Chenab (J-C) system. The mean canal water requirements and associated share allocations that need to be met from stored water are shown in Figure 2c, highlighting the typical quantity and timing differences between the two systems that must be balanced. Figure 2d shows the cumulative stored water requirements to meet actual canal deliveries and share requirements throughout a year, in addition to current available storage volume in Tarbela and Mangla. The figure shows that the IBIS typically needs water stored in both Tarbela and Mangla reservoirs to meet the Rabi seasonal demand. Without considering delivery losses, the current available storage is only sufficient to meet typical Rabi season demand patterns, provided both reservoirs are operated in harmony with an extensive link canal system to capture and distribute water equitably. The current available storage volume is insufficient to meet the water sharing requirements of the 1991 Water Apportionment Accord (the Accord).
A complex seasonal inflow forecast, and storage operational planning process is undertaken on the IBIS at the start of every season to meet Indus and J-C system allocation demand patterns, while maintaining equitable sharing between provinces and then command canals as laid out in the Accord. The Accord is the agreement that describes the distribution of the water resources of the Indus between the 4 provinces of Punjab, Sindh, Balochistan and Khyber Pakhtunkhwa (KP). In 14 paragraphs, the Accord outlines how the water resources are to be shared on a seasonal and 10-day basis [38]. There are 3 paragraphs that are relevant to the sharing arrangements (Appendix A)-Paragraphs 14(b), 2 and 4. Paragraph 14(b) provides seasonal totals based on observed average use for the period 1977-1982. Paragraph 2 considers accepted water distributional principles and details shares on a provincial and seasonal basis (Figure 2a,b) and Paragraph 4 is used to balance supplies above Paragraph 2 from floods and future storages. The precise interpretation and subsequent practical implementation of the Accord paragraphs are not documented [20,21] and extensive consultation with Indus River System Authority (IRSA) operators to elicit these operational rules was required for this modelling study. This component is essential for any simulation of allocation and subsequent distribution of water within the IBIS and is outlined below.
Water 2021, 13, x FOR PEER REVIEW 5 of 21 erators to elicit these operational rules was required for this modelling study. This component is essential for any simulation of allocation and subsequent distribution of water within the IBIS and is outlined below. The allocation process is undertaken at the start of Kharif and Rabi seasons, with the Kharif season being sub-divided into the early (1st April-10th June) and late Kharif (11th June-30th September). The primary objective of the allocation process is to maximise resource allocation while maintaining 10-day provincial shares, as specified in the Accord. This is done while balancing water between two supply sub-systems-Indus and J-C. Operationally, the provincial shares determined at the start of the water year are managed during the season with adjustments, as the observed inflows differ from forecasts. The aim of the adjustments is to achieve the total seasonal provincial share proportions calculated at the start of the season. Provincial sharing to command canals is based on defined ratios that vary due to operational constraints, such as canal closures and is undertaken at a provincial level.
There are 7 components to the resource assessment and allocation process ( Figure 3); (i) 10-day inflow forecasts at 5 rim stations (Kabul, Indus, Jhelum, Chenab and Eastern Rivers); (ii) unaccounted loss/gain estimation for J-C and Indus systems; (iii) minimum and maximum 10-day storage level targets for Tarbela and Mangla reservoirs; (iv) reservoir balancing; (v) resource balancing across 10-day periods and supply systems; (vi) provincial resource sharing and (vii) sharing to provincial command canals.
Currently, IRSA uses a probabilistic method for the inflow forecast in the western rivers (Kabul, Indus, Jhelum and Chenab). The eastern river forecasts (Ravi and Sutlej) are based on recent observed inflow history. Other forecast methods are being trialed within the IBIS and provision has been made in the model to include these when they are adopted. The allocation process is undertaken at the start of Kharif and Rabi seasons, with the Kharif season being sub-divided into the early (1 April-10 June) and late Kharif (11 June-30 September). The primary objective of the allocation process is to maximise resource allocation while maintaining 10-day provincial shares, as specified in the Accord. This is done while balancing water between two supply sub-systems-Indus and J-C. Operationally, the provincial shares determined at the start of the water year are managed during the season with adjustments, as the observed inflows differ from forecasts. The aim of the adjustments is to achieve the total seasonal provincial share proportions calculated at the start of the season. Provincial sharing to command canals is based on defined ratios that vary due to operational constraints, such as canal closures and is undertaken at a provincial level.
There are 7 components to the resource assessment and allocation process ( Figure 3); (i) 10-day inflow forecasts at 5 rim stations (Kabul, Indus, Jhelum, Chenab and Eastern Rivers); (ii) unaccounted loss/gain estimation for J-C and Indus systems; (iii) minimum and maximum 10-day storage level targets for Tarbela and Mangla reservoirs; (iv) reservoir balancing; (v) resource balancing across 10-day periods and supply systems; (vi) provincial resource sharing and (vii) sharing to provincial command canals.
Currently, IRSA uses a probabilistic method for the inflow forecast in the western rivers (Kabul, Indus, Jhelum and Chenab). The eastern river forecasts (Ravi and Sutlej) are based on recent observed inflow history. Other forecast methods are being trialed within the IBIS and provision has been made in the model to include these when they are adopted. Water 2021, 13, x FOR PEER REVIEW 7 of 21

Modelling Platform
Considering the physical complexities and rule-based water distribution procedures in IBIS, the eWater Source framework was adopted to construct IRSM, as it contains mechanisms to simulate planning and decision-making aspects of seasonal water demand. The eWater Source modelling framework allows catchment-scale rainfall runoff modelling (using a geographic interface) and rule-based mass balance modelling through a schematic node-link network with flow routing [39][40][41][42]. The eWater Source modelling framework is used extensively throughout Australia [43], as well as in major international river basins. For this study, observed rim flows were used, rather than simulated upper Indus flows. For the lower Indus basin, a node-link structure was constructed to facilitate simulation of the water flow through the key IBIS infrastructure, whereas a series of custom functions was implemented to simulate the water-sharing system. As shown in Figure 4, there are 52 main canal commands in the IBIS with offtakes from barrages supplying 72 irrigation areas. The conceptualisation of the supply-based system includes the physical dimensions of the barrage and canal infrastructure that is needed to constrain water deliveries. Full details of barrage and offtake canal configurations can be found in Stewart et al. [44]. Groundwater can also be drawn upon to meet irrigation shortfalls if the groundwater is of sufficient quality. The IRSM is described in detail in Stewart et al. [44]. For this study, a more accurate representation of the seasonal allocation planning procedure of IRSA was required that accounted for the change in available storage in the planning phase of seasonal operations and allocation; therefore, an external seasonal forecast as described in Section 2.1.4 was generated to drive the provincial water-sharing functions in the IRSM, which then applied these forecasts and share ratios to canal commands in the model domain to drive the water ordering system functions in the IRSM.
The availability of data described in Appendix B influenced the conceptualisation of the model in terms of the storages, barrages and reaches included in the model. Demand The western river forecasts are based on using the respective previous season inflow volume as an indicator of the forthcoming season inflows. Historic (from 1976) previous seasons' inflow volumes that are within 5% of the current previous season inflow provide a list of forthcoming season probabilities that are averaged to the nearest 5 percentile to provide a forthcoming season probability (Rabi, early and Late Kharif). Then, using probabilities of 10% more and less, the 10-day inflow values are selected for the respective minimum and maximum inflow forecasts.
The eastern river forecast is based on the last five years of corresponding seasonal inflow, with the largest and smallest being the respective maximum and minimum 10-day inflow sequence. The eastern river component of the overall forecast is small and is dependent on river operations outside of Pakistan. Although other forecasting methods may be available, the last-five-years method has been adopted as suitable by IRSA and has therefore been implemented in this modelling study.
Unaccounted loss/gain proportions for Rabi, early and late Kharif are worked out separately for the Indus and J-C systems and are based on the average over the most recent 10 years. For the Indus system, it is the change in Tarbela and Chashma storage volume plus Indus canal withdrawals and below Kotri flows, minus the total Indus inflow (Kabul and Indus inflows plus river releases from Rasul, Qadirabad, Balloki and Sulemanki barrages), divided by the Indus system inflow. The J-C system is the change in Mangla storage volume plus J-C canal withdrawals and river releases from Rasul, Qadirabad, Balloki and Sulemanki barrages minus the J-C inflow (Jhelum, Chenab, Ravi and Sutlej inflows), divided by the J-C system inflow.
There are seasonal 10-day minimum and maximum storage level targets for Tarbela and Mangla. The Kharif filling targets ensure a minimum Tarbela filling rate (between 0 and 0.03048 m/day) and a maximum filling rate of 0.3048 m/day past 460 m reservoir elevation. Mangla is constrained to release at least 227 m 3 /s for Rasul hydro canal operation, with an 80% filling target by the 30 June. Tarbela is to be filled by the 20th and 31st of August for respective maximum and minimum inflow scenarios and Mangla is to be filled by the end of August. There is a 10% drawdown of Tarbela and Mangla volume at the end of September. Rabi has storage drawdown targets for the end of December (10%) and March 20th (empty). Filling rates and storage targets consider reservoir capacity reductions through time due to sedimentation.
Provincial sharing is based on 10-day patterns specified in the 1991 Accord paragraphs. Paragraph 14(b) patterns are used in Rabi and Paragraph 2 patterns are used in Kharif. The calculation of the available resource is done for each of the systems. Any shortages to meet these constraints are equalised between sub-systems and across 10-day periods by managing outflows from the J-C to Indus system. The method ensures that at the end of the allocation process, flows below Kotri Barrage (the end of the system) are minimised, reservoirs have identical volume proportions for each 10-day period and Indus and J-C systems have identical shortages for users in these sub-systems.
Once the 10-day available resource is determined, the water is shared firstly to KP and Balochistan, with the remaining resource shared to Punjab and Sindh according to a 3 tier approach on a 10-day basis and reconciled for the season. If the resource is less than the Paragraph 14(b) Punjab and Sindh total share, then it is shared according to paragraph 14(b); if it is between the Paragraph 14(b) and 2 Punjab and Sindh total share, the resource above 14(b) is shared based on an interpolated ratio until it reaches the Paragraph 2 value, whereby it is shared according to Paragraph 2. The provincial shares are then distributed to canals based on provincial canal shares. When surplus water is available above Paragraph 2, provinces can request additional supplies according to Paragraph 4.

Model Conceptualisation 2.2.1. Modelling Platform
Considering the physical complexities and rule-based water distribution procedures in IBIS, the eWater Source framework was adopted to construct IRSM, as it contains mechanisms to simulate planning and decision-making aspects of seasonal water demand. The eWater Source modelling framework allows catchment-scale rainfall runoff modelling (using a geographic interface) and rule-based mass balance modelling through a schematic node-link network with flow routing [39][40][41][42]. The eWater Source modelling framework is used extensively throughout Australia [43], as well as in major international river basins. For this study, observed rim flows were used, rather than simulated upper Indus flows. For the lower Indus basin, a node-link structure was constructed to facilitate simulation of the water flow through the key IBIS infrastructure, whereas a series of custom functions was implemented to simulate the water-sharing system. As shown in Figure 4, there are 52 main canal commands in the IBIS with offtakes from barrages supplying 72 irrigation areas. The conceptualisation of the supply-based system includes the physical dimensions of the barrage and canal infrastructure that is needed to constrain water deliveries. Full details of barrage and offtake canal configurations can be found in Stewart et al. [44]. Groundwater can also be drawn upon to meet irrigation shortfalls if the groundwater is of sufficient quality. The IRSM is described in detail in Stewart et al. [44]. For this study, a more accurate representation of the seasonal allocation planning procedure of IRSA was required that accounted for the change in available storage in the planning phase of seasonal operations and allocation; therefore, an external seasonal forecast as described in Section 2.1.4 was generated to drive the provincial water-sharing functions in the IRSM, which then applied these forecasts and share ratios to canal commands in the model domain to drive the water ordering system functions in the IRSM.
The availability of data described in Appendix B influenced the conceptualisation of the model in terms of the storages, barrages and reaches included in the model. Demand data constrained the model to the command canal level, with groundwater data identifying areas of conjunctive use.

Water Allocation and Distribution System
The Accord water-sharing is implemented in Source via a series of functions [45] that take input from the provincial scale forecasting and sharing, as described in Section 2.1.4, which is currently undertaken outside of the Source modelling framework in the Pakistan Water Apportionment Accord (WAA) tool [46,47]. Previous IRSM simulations [44] implemented a simplified provincial allocation system entirely in Source functions; however, this did not capture the iterative storage balancing required or the updating of the probabilistic forecasting approach as the model progresses through each year of simulation.
In the model described in this paper, the provincial forecasting and sharing is undertaken outside of Source at the beginning of each season (Kharif and Rabi). Rim station inflows from each subsequent year of simulation are used to rebuild the probability tables underpinning the seasonal forecasting and generate forecasts for the forthcoming season. The resultant seasonal provincial allocation is then input to Source, where a series of func-

Water Allocation and Distribution System
The Accord water-sharing is implemented in Source via a series of functions [45] that take input from the provincial scale forecasting and sharing, as described in Section 2.1.4, which is currently undertaken outside of the Source modelling framework in the Pakistan Water Apportionment Accord (WAA) tool [46,47]. Previous IRSM simulations [44] implemented a simplified provincial allocation system entirely in Source functions; however, this did not capture the iterative storage balancing required or the updating of the probabilistic forecasting approach as the model progresses through each year of simulation.
In the model described in this paper, the provincial forecasting and sharing is undertaken outside of Source at the beginning of each season (Kharif and Rabi). Rim station inflows from each subsequent year of simulation are used to rebuild the probability tables underpinning the seasonal forecasting and generate forecasts for the forthcoming season. The resultant seasonal provincial allocation is then input to Source, where a series of functions disaggregate the forecast allocations to create time-series demands at individual canal commands, implemented at minimum flow requirement nodes (Figure 4 inset). The Source ordering priority pathways ( Figure 4) are then used to distribute the collated provincial allocations to reservoirs, adding travel time and associated instream losses along the way to ensure that sufficient water can be released. If insufficient water is available to be released at any given time to meet orders, all orders are proportionally reduced. The link canals of C-J Link and T-P link offer the choice of switching ordering priority pathways between Tarbela and Mangla reservoirs. Functions based on average behavior are installed to look up relative reservoir volumes and select the pathway to attempt to harmonise the storages. Finally, reservoir, barrage and offtake canal and link canal capacities provide a further limiting capacity for both the ordering network and the delivery network.

Calibration and Validation Methodology
A stepwise manual calibration approach was used for the calibration of individual reaches to account for river routing and unaccounted losses and gains. For operational decisions and water allocation and sharing aspects, interview/data assimilation techniques were adopted [48]. The identification of priority flow pathways ( Figure 4) combined with water allocation to canal commands has by far the largest influence on the performance of this model at barrages and canal command offtakes. Although automated calibration techniques could potentially have been used to determine these priority pathways and water-sharing arrangements, interviews with the appropriate water managers was considered a more representative method of replicating actual behavior. Therefore, the calibration approach was deemed to be appropriate for the operational decision-making.
The IRSM was calibrated using streamflow and level-volume data from 2007 to 2013. This period represented approximately half of the detailed data set available to the study team, including daily canal water delivery data. The calibration period contained both wet and dry years and was considered long enough to identify the priority flow pathways enacted by river managers. Model validation was for the period of 2002 to 2007 and an extended simulation period of 1990 to 2013 was used for long-term model validation.
Calibration metrics used for flows at barrages and canal command heads were the daily Nash-Sutcliffe coefficient of efficiency (NSE) [49] and overall bias (%). For river routing, piecewise linear routing [50,51] is used. For each reach with observed upstream and downstream flows a piecewise flow-versus-travel-time relationship is assigned based on observed peak travel times. Following assignment of the piecewise table, reach unaccounted differences were estimated by forcing upstream flows in the reach and routing these flows to the downstream observation location and comparing the observed and modelled flow duration curves. The mean difference between the curves was used to form a flow-versus-unaccounted-difference relationship for the reach. Reach routing and unaccounted differences were then combined and optimised for overall bias and daily NSE correlation. Reach models were then combined (stitched) in the IRSM and fine-tuned where appropriate. All results in this paper reflect stitched model results.
For operational decisions, interviews with IRSA, Water and Power Development Authority (WAPDA), Punjab and Sindh Irrigation Department personnel were consulted to determine preferred flow pathways and storage/barrage operational data [52]. Water orders generated from the seasonal allocation process are passed from most canal commands up the river/canal pathways to Chashma and Mangla. These orders accumulated travel times and losses so that the correct amount of water was released to arrive at the destination in time. For Eastern Canal Commands, orders can only be met by Marala (Figure 4). The identified water order priority pathways give the model options to source water from different storages in the model, depending on availability. Orders passed up through T-P link and C-J link canals are dependent on the modelled available storage levels in Mangla and Tarbela, with the order directed to the storage with the most available water. Water arriving from unregulated tributaries or canal branches are used to reduce orders to storages.
Depending on inflows, barrage water levels are generally allowed to fluctuate in the model throughout the season; however, some barrages have specific operational heads assigned throughout the year. In such cases, the operational head is set based on the analysis of historic water levels at the barrage for different times of the year.
For water allocation and sharing, the processes adopted by IRSA, as described in Section 2.1.4, were adopted in the model and comparisons between overall modelled and observed share ratios were used to assess model performance of water allocation and sharing components, in addition to overall bias between observed and modelled canal extractions.

Sedimentation Scenario
The observed reduction in storage volume resulting from sedimentation through time for Tarbela and Mangla is shown in Figure 5. To investigate these impacts of reservoir sedimentation on water resource availability, two scenarios are considered. The 2020 scenario represents current infrastructure and water-sharing rules with current levels of sedimentation, Tarbela and Mangla capacities are respectively 7.38 km 3 (in 2020) and 9.08 km 3 (2018). The 1990 scenario considers the current conditions of the model with 1990 levels of sedimentation.
water. Water arriving from unregulated tributaries or canal branches are used to reduce orders to storages.
Depending on inflows, barrage water levels are generally allowed to fluctuate in the model throughout the season; however, some barrages have specific operational heads assigned throughout the year. In such cases, the operational head is set based on the analysis of historic water levels at the barrage for different times of the year.
For water allocation and sharing, the processes adopted by IRSA, as described in Section 2.1.4, were adopted in the model and comparisons between overall modelled and observed share ratios were used to assess model performance of water allocation and sharing components, in addition to overall bias between observed and modelled canal extractions.

Sedimentation Scenario
The observed reduction in storage volume resulting from sedimentation through time for Tarbela and Mangla is shown in Figure 5. To investigate these impacts of reservoir sedimentation on water resource availability, two scenarios are considered. The 2020 scenario represents current infrastructure and water-sharing rules with current levels of sedimentation, Tarbela and Mangla capacities are respectively 7.38 km 3 (in 2020) and 9.08 km 3 (2018). The 1990 scenario considers the current conditions of the model with 1990 levels of sedimentation.

Scenario 1: Current Conditions (2020)
This scenario adopts provincial shares under 2020 reservoir operational capacities with the current sedimentation volume using the historic inflow regime from 1990-2017. Using this configuration of sedimentation, the Tarbela and Mangla capacities are respectively 7.38 (in 2020) and 9.08 km 3 (2018).

Scenario 2: Current Conditions with 1990 Sediment Levels
This scenario adopts provincial shares under 2020 water reservoir capacities but with 1990 sedimentation volumes (i.e., greater operational capacity of the reservoirs) using the historic inflow regime from 1990-2017. The capacities of Tarbela and Mangla for this scenario are 9.50 km 3 and 9.36 km 3 , respectively. Note that this scenario does not represent 1990 levels of infrastructure, e.g., Mangla capacity prior to enlargement in 2009.

Scenario 1: Current Conditions (2020)
This scenario adopts provincial shares under 2020 reservoir operational capacities with the current sedimentation volume using the historic inflow regime from 1990-2017. Using this configuration of sedimentation, the Tarbela and Mangla capacities are respectively 7.38 (in 2020) and 9.08 km 3 (2018).

Scenario 2: Current Conditions with 1990 Sediment Levels
This scenario adopts provincial shares under 2020 water reservoir capacities but with 1990 sedimentation volumes (i.e., greater operational capacity of the reservoirs) using the historic inflow regime from 1990-2017. The capacities of Tarbela and Mangla for this scenario are 9.50 km 3 and 9.36 km 3 , respectively. Note that this scenario does not represent 1990 levels of infrastructure, e.g., Mangla capacity prior to enlargement in 2009. Figure 6 shows the mean annual unaccounted difference for river reaches, where both routing and unaccounted difference were modelled, and highlights the lower Indus reaches where the largest unaccounted differences were. These reaches, however, are generally not those with the highest proportional unaccounted differences, which are generally associated with the eastern rivers. Figure 6 shows the mean annual unaccounted difference for river reaches, where both routing and unaccounted difference were modelled, and highlights the lower Indus reaches where the largest unaccounted differences were. These reaches, however, are generally not those with the highest proportional unaccounted differences, which are generally associated with the eastern rivers. In matching observed flows at barrages, daily operational decisions and the setting of preferential flow paths for water orders were found to have the greatest influence on model performance. After accounting for these operational decisions in the model based on average or most common behaviour, the mean annual relative error (bias)-expressed as a percentage of mean annual river flow for each reach (denoted by downstream barrage location)-is presented in Figure 7 and associated daily NSE, shown in Figure 8. In matching observed flows at barrages, daily operational decisions and the setting of preferential flow paths for water orders were found to have the greatest influence on model performance. After accounting for these operational decisions in the model based on average or most common behaviour, the mean annual relative error (bias)-expressed as a percentage of mean annual river flow for each reach (denoted by downstream barrage location)-is presented in Figure 7 and associated daily NSE, shown in Figure 8.

Physical System Performance
The river section results showed that: • River reaches with the highest flows had the lowest relative error percentage and the best daily NSE correlation.

•
The daily NSE correlation was generally high but diminished with increasing distance from the original water source (i.e., Sidhnai and Islam Barrages). The poor performance at Rasul Barrage was due to a poor correlation with observed releases from Mangla Reservoir.
• Barrages with low inflows and most of the water being diverted to canals such as Sidhnai, Balloki, Khanki and Punjnad performed the worst because small or moderate errors in either the daily correlation associated with the upstream barrage or volume were amplified in the small downstream releases. • Simulated flows exhibited a positive bias upstream of Guddu barrage at Punjnad and Sidhnai. This performance issue was associated with the distribution of water between Indus and J-C systems via TP link canal and corresponding positive bias at Balloki. A greater proportion of water was distributed through eastern link canals in the model, resulting in a negative bias at Guddu Barrage, which was exacerbated further downstream.

•
The simulated Indus flows exhibited an additional negative bias downstream of Sukkur Barrage, which may be due to confusion around whether Balochistan deliveries and the Karachi water supplies from Kotri were included in observed Sindh withdrawals.  The river section results showed that: • River reaches with the highest flows had the lowest relative error percentage and the best daily NSE correlation.

•
The daily NSE correlation was generally high but diminished with increasing distance from the original water source (i.e., Sidhnai and Islam Barrages). The poor performance at Rasul Barrage was due to a poor correlation with observed releases from Mangla Reservoir.

•
Barrages with low inflows and most of the water being diverted to canals such as Sidhnai, Balloki, Khanki and Punjnad performed the worst because small or moderate errors in either the daily correlation associated with the upstream barrage or volume were amplified in the small downstream releases. • Simulated flows exhibited a positive bias upstream of Guddu barrage at Punjnad and Sidhnai. This performance issue was associated with the distribution of water be-  The river section results showed that: • River reaches with the highest flows had the lowest relative error percentage and the best daily NSE correlation.

•
The daily NSE correlation was generally high but diminished with increasing distance from the original water source (i.e., Sidhnai and Islam Barrages). The poor performance at Rasul Barrage was due to a poor correlation with observed releases from Mangla Reservoir.

•
Barrages with low inflows and most of the water being diverted to canals such as Sidhnai, Balloki, Khanki and Punjnad performed the worst because small or moderate errors in either the daily correlation associated with the upstream barrage or volume were amplified in the small downstream releases. • Simulated flows exhibited a positive bias upstream of Guddu barrage at Punjnad and Sidhnai. This performance issue was associated with the distribution of water be-  Figure 9 shows that there was an average canal error of 5%, with the worst performance at Taunsa T-P Link at 15%. This is the canal that is used to balance usage between Mangla and Tarbela and is difficult to match daily operational decisions. The data shows that most canals exhibited less than 10% bias, which is well within the measurement error. Figure 10 shows that the maximum percentage volumetric errors in Tarbela and Mangla were as large as 50%. The figure also shows that an error in one reservoir is typically offset by a similar magnitude negative error in the other. This finding indicates that the model balances the reservoirs differently compared with observed operational behaviour. When combined, the overall storage volume error oscillates between positive (maximum 31%) and negative (minimum −30%), indicating no bias toward retaining or releasing too much water. As the storage capacities of the dams are much less than the mean annual flow in the rivers, the dam storage volumetric error, expressed as a percentage of mean annual flow is relatively small. For the calibration period, the Mangla, Tarbela and combined respective storage volume errors relative to inflows are −6.4% to 8.8%, −5.2% to 3.6% and −4.3% to 4.5% for the combined system.
Water 2021, 13, x FOR PEER REVIEW 13 of 21 Figure 9 shows that there was an average canal error of 5%, with the worst performance at Taunsa T-P Link at 15%. This is the canal that is used to balance usage between Mangla and Tarbela and is difficult to match daily operational decisions. The data shows that most canals exhibited less than 10% bias, which is well within the measurement error.  Figure 10 shows that the maximum percentage volumetric errors in Tarbela and Mangla were as large as 50%. The figure also shows that an error in one reservoir is typically offset by a similar magnitude negative error in the other. This finding indicates that the model balances the reservoirs differently compared with observed operational behaviour. When combined, the overall storage volume error oscillates between positive (maximum 31%) and negative (minimum −30%), indicating no bias toward retaining or releasing too much water. As the storage capacities of the dams are much less than the mean annual flow in the rivers, the dam storage volumetric error, expressed as a percentage of mean annual flow is relatively small. For the calibration period, the Mangla, Tarbela and combined respective storage volume errors relative to inflows are −6.4% to 8.8%, −5.2 to 3.6% and −4.3% to 4.5% for the combined system. The model is configured to match average behaviour and does not capture all daily operational decisions and consequently an exact match with observed behaviour is not possible. Most of the model's uncertainty is caused by daily operational decisions. The model is unable to simulate subjective decisions, but it tried to replicate the average of balancing of reservoirs and distributions of water between systems across the historical period. As the inflows are considerably larger than the stored resource, the allocation process is more sensitive to the inflows than where the stored resource resides. Therefore, it that most canals exhibited less than 10% bias, which is well within the measurement error.  Figure 10 shows that the maximum percentage volumetric errors in Tarbela and Mangla were as large as 50%. The figure also shows that an error in one reservoir is typically offset by a similar magnitude negative error in the other. This finding indicates that the model balances the reservoirs differently compared with observed operational behaviour. When combined, the overall storage volume error oscillates between positive (maximum 31%) and negative (minimum −30%), indicating no bias toward retaining or releasing too much water. As the storage capacities of the dams are much less than the mean annual flow in the rivers, the dam storage volumetric error, expressed as a percentage of mean annual flow is relatively small. For the calibration period, the Mangla, Tarbela and combined respective storage volume errors relative to inflows are −6.4% to 8.8%, −5.2 to 3.6% and −4.3% to 4.5% for the combined system. The model is configured to match average behaviour and does not capture all daily operational decisions and consequently an exact match with observed behaviour is not possible. Most of the model's uncertainty is caused by daily operational decisions. The model is unable to simulate subjective decisions, but it tried to replicate the average of balancing of reservoirs and distributions of water between systems across the historical period. As the inflows are considerably larger than the stored resource, the allocation process is more sensitive to the inflows than where the stored resource resides. Therefore, it The model is configured to match average behaviour and does not capture all daily operational decisions and consequently an exact match with observed behaviour is not possible. Most of the model's uncertainty is caused by daily operational decisions. The model is unable to simulate subjective decisions, but it tried to replicate the average of balancing of reservoirs and distributions of water between systems across the historical period. As the inflows are considerably larger than the stored resource, the allocation process is more sensitive to the inflows than where the stored resource resides. Therefore, it is not sensitive to average operational behaviour assumptions. However further refinement of storage target levels and order distribution between Indus and J-C systems could improve modelled storage performance in the model, but this would be outside of the agreed average operational behaviour. On this basis, the average behaviour over the modelling period was reasonably good and is a considerable improvement on existing models that do not accurately represent the sharing and distribution rules.

Water Allocation and Sharing Performance
The modelled provincial flow proportions are a close match to historic allocations across both the calibration and validation period. Figure 11 shows that the model replicated the average water-sharing between Sindh and Punjab within 0.1% with seasonal error (not shown) less than 1%. that do not accurately represent the sharing and distribution rules.

Water Allocation and Sharing Performance
The modelled provincial flow proportions are a close match to historic allocations across both the calibration and validation period. Figure 11 shows that the model replicated the average water-sharing between Sindh and Punjab within 0.1% with seasonal error (not shown) less than 1%.     Figure 12 shows the time series of modelled and observed total canal withdrawals for Punjab and Sindh provinces. The match is reasonably good, considering that the model does not consider specific operational decisions such as temporary canal closures for maintenance or closures during flood periods (i.e., 2010).

Water Allocation and Sharing Performance
The modelled provincial flow proportions are a close match to historic allocations across both the calibration and validation period. Figure 11 shows that the model replicated the average water-sharing between Sindh and Punjab within 0.1% with seasonal error (not shown) less than 1%.

Reservoir Sedimentation Impacts
The operational capacities of Tarbela and Mangla reservoirs are being reduced due to sedimentation; however, the water-sharing policies are based on fixed shares. The IRSM model was run with current (2020) level-volume relationships for Tarbela and Mangla and with previous relationships for 1990 and run for the historic period of 1990-2017 (26 complete years) to assess the change in available resources and the shares to Punjab and Sindh. Table 1 shows the modelled change in seasonal shares and resources for the two scenarios. Balochistan and KP have high security shares that do not change. Table 1 shows that for a 2.4 km 3 reduction in storage capacity there is an annual reduction of 0.95 km 3 of supplied resource, which is predominantly in Rabi. In early Kharif there is a slight increase in resources for Punjab and Sindh in the 2020 scenario. This is because less water is being set aside for reservoir filling and can therefore be allocated.
However, in late Kharif, where spilling is increased due to sedimentation, Sindh may be able to opportunistically take unregulated water to fulfil share requirements. In percentage terms, the Rabi supply for Punjab and Sindh combined is respectively reduced on average by between 4% and 5% across equivalent years of simulation ( Figure 13), whereas the average increase in Early Kharif and Kharif flows is less than 1% across equivalent years. In some years, the modelling indicates larger Kharif percentage changes in canal deliveries, but generally not lower canal deliveries. All 26 years of model simulation show a decrease in total Rabi deliveries, with the largest percentage reductions in deliveries occurring in the late October to early December period.
with previous relationships for 1990 and run for the historic period of 1990-2017 (26 complete years) to assess the change in available resources and the shares to Punjab and Sindh. Table 1 shows the modelled change in seasonal shares and resources for the two scenarios. Balochistan and KP have high security shares that do not change.  Table 1 shows that for a 2.4 km 3 reduction in storage capacity there is an annual reduction of 0.95 km 3 of supplied resource, which is predominantly in Rabi. In early Kharif there is a slight increase in resources for Punjab and Sindh in the 2020 scenario. This is because less water is being set aside for reservoir filling and can therefore be allocated. However, in late Kharif, where spilling is increased due to sedimentation, Sindh may be able to opportunistically take unregulated water to fulfil share requirements.
In percentage terms, the Rabi supply for Punjab and Sindh combined is respectively reduced on average by between 4% and 5% across equivalent years of simulation ( Figure  13), whereas the average increase in Early Kharif and Kharif flows is less than 1% across equivalent years. In some years, the modelling indicates larger Kharif percentage changes in canal deliveries, but generally not lower canal deliveries. All 26 years of model simulation show a decrease in total Rabi deliveries, with the largest percentage reductions in deliveries occurring in the late October to early December period.

Discussion
The modelled scenario shows that the reduction in storage capacity does not equate to an equivalent reduction in canal deliveries. The Rabi season canal deliveries see the greatest percentage decline between the 1990 and 2020 scenarios of between 4% and 5% or approximately 1.26 km 3 of canal deliveries per year. This difference is due to the system harvesting opportunistic (unregulated) water to meet supply share requirements. The most heavily impacted period of the Rabi is late October to early December, i.e., the crop sowing period, when stored water supplies run out. The month of January, however, coincides with canal closures and is associated with the lowest canal delivery period. During February and March, opportunistic water can be taken (provided the water is not stored), further reducing the difference between scenarios at this time of the year. Although many studies may look at potential changes to inflows to the IBIS from climate change or potential upstream development, which are uncertain, this study has focused on a historic decrease in available storage through sedimentation, which has shown remarkably consistent trends that will continue. This would indicate that without further augmentation of system storage, the Rabi season supplies will continue to be further impacted in the future. An essential component of assessing these future supply scenarios will be the potential changes to inflows and associated sediment load because of climate change. This study is limited in that it does not yet consider these potential future impacts. The model also does not attempt to identify the associated increase in crop failures or critical supply thresholds; however, these limitations may be considered in future modelling studies.
The unaccounted flow differences highlighted in Figure 6 are a large component of the water balance and are consequently an area of concern for water managers. The model calibration uses average unaccounted difference estimates for each of the reaches and calibration for overall bias in a reach did not identify a need for changes to the nonlinear flow vs. unaccounted difference tables over time. However, considerable debate is ongoing in Pakistan over possible trends in unaccounted differences. We suggest that perceived changes to unaccounted differences may in part be explained by changes in internal and external flows which are captured by the model through the non-linear flow vs. unaccounted difference tables. In many reaches of the IRSM, the derived unaccounted difference term diminishes as a proportion of flow as flows increase until out-of-bank (flood) flow conditions are met. These relationships suggest that in low flow years, unaccounted differences as a proportion of flow rate may increase, although further research is required to confirm this.
Few river modelling frameworks have sufficient flexibility to include a seasonal planning phase as complex as that used in the IBIS. The IRSM originally contained functions that approximated the seasonal planning calculations; however, to capture a more complete process, the seasonal forecast calculations were undertaken outside of the IRSM. This is a key limitation of the current version of IRSM and work is currently underway to better integrate the planning phase of the IBIS through a software plugin to Source. The current model has, however, demonstrated key operational aspects of the IBIS that are not addressed by other models, namely, seasonal water allocation and distribution to canal commands and river operations to link canals and canal commands. The IRSM captures previously undocumented operational rules and expert knowledge essential to the representation of the IBIS.
Finally, other current limitations of the IRSM include no allowance for urban population growth over time and incomplete representation of hydropower production. In many cases, urban surface water supplies are small in relation to unaccounted differences and this does not impact model accuracy. However, Karachi extracts a significant amount of water from the Indus and has significant population growth [53]. Future model studies could be improved by a better conceptualisation of the Karachi water supply. Hydropower production has been explored by the model; however, as Figure 10 demonstrates, hydropower production may prove difficult to match without first achieving a closer match with reservoir storage operations.

Conclusions
In response to the Water Sector Task Force Report of 2012, the IRSM has been built to replicate both the physical and management (allocation, sharing and distribution) processes in the Pakistan IBIS. The purpose of this study was to assess the capability and suitability of the IRSM to model the system's operation and assess the impacts of reservoir sedimentation on water resource security and sharing.
Model calibration results for the period 2002-2012 show that provincial allocation ratios agree with historic water deliveries (0.1% error in sharing) and that water delivered to barrages and then on to canal commands is also in agreement with recorded data (major barrages with a daily NSE of 0.7). Storage behaviour simulated by the model indicates some scope for improving the balancing of stored water resources between storages. The average maximum volume errors with respect to the mean annual inflow for Tarbela and Mangla are respectively 5.2% and 8.8%. This performance demonstrates that the IRSM is capable of modelling water-sharing and allocations to replicate the historic allocation of water in the IBIS. Furthermore, the model can also be used to investigate the unaccounted water in the system and to ensure necessary efficiency in the water distribution. Based on this performance, the model is suitable to investigate the impact of climate change and reservoir sedimentation on canal command water balance analysis and future provincial water security in Pakistan [53,54]. In addition to impact assessment, this modelling framework can be enhanced to investigate the whole of the system and canal-command-scale current and future water planning and management decisions.
The reservoir sedimentation study showed that a 2.3 km 3 reduction in storage volume equates to a 1.26 km 3 reduction in water supplied in the Rabi season. It also showed that the impact is not equitable between provinces, with KP and Balochistan not being affected in similar way to Punjab and Sindh. This problem will exacerbate over time as the storages continue to fill with sediment and changes to water-sharing arrangements will be required in response to this. This study found that the potential reduction in Rabi canal deliveries resulting from reservoir sedimentation were partly offset by small increases in Kharif canal deliveries; however, the model also indicated that some parts of the Rabi season experienced more severe shortages that may reach critical points in the future as reservoir sedimentation continues.

Data Availability Statement:
The data that support the findings of this study are available from the relevant Pakistan Government organisation. Restrictions apply to the availability of these data, which were used under license for this study. Data are available from the authors with the permission of relevant Pakistan Government organisation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, and in the collection, analyses, or interpretation of data. * Adopted from Government of Pakistan [38], and Khan [55].