Early twenty-first century satellite-driven irrigation performance in the world’s largest system: Pakistan’s Indus Basin irrigated system

Evaluating irrigation performance in large systems is often limited by the availability of reliable water use data. Satellite-driven actual evapotranspiration (ETa) estimates are used herein as water use surrogates to assess the year-to-year inter-seasonal irrigation performance in 46 canal commands in the Indus Basin irrigated system (IBIS), the largest in the world (∼160 000 km2). The accuracy and reliability of the ETa estimates are verified using two previously published locally adjusted satellite-driven ETa estimates, as well as field ETa estimates. Inter-seasonal variability (canal command water use in time) and equity (inter- and intra-canal command water use) are assessed from 2000 to 2018 using violin-plots time-series for the two irrigation seasons, the wet ‘Kharif’ and dry ‘Rabi’. The violin-plots probability density functions are used to assess intra-canal command equity; and their seasonal time-series to assess inter-seasonal variability. The long-term multi-year assessment conducted here, the first for the IBIS using consistent satellite-driven ETa time-series, shows that canal commands with ready access to groundwater exhibit more equity and less inter-seasonal variability when compared to canal commands chiefly reliant on surface water supplies; with the latter showing intra-canal command inequities between head-end and tail-end irrigated areas. Also, ETa in canal commands is mostly slightly increasing and there is low inter-seasonal variability in both irrigation seasons, except for two canal command at the system-end, which show higher inter-seasonal variability and inequity than their upstream counterparts. The methods employed here can be used in large irrigated systems elsewhere to assess ongoing irrigation performance and to verify results of targeted (non)structural irrigation management.


Introduction
Increased reliance of surface and groundwater to meet the demand for food and other agricultural products is placing water resources on an unsustainable trajectory in many parts of the world (Gleeson and Wada 2013, Wada andBierkens 2014). The largest irrigated system in the world, the Indus Basin irrigated system (IBIS, ∼160 000 km 2 , FAO 2012) in Pakistan is arguably on this trajectory; with declining groundwater levels and salinity impacting water availability and crop productivity (Briscoe and Qamar 2005, Khan et al 2017, Kirby et al 2017. The IBIS is the lifeline of Pakistan, serving ∼90 000 farmers and providing food and economic security for ∼207 million people (FAO 2012). The Chenab, Jhelum, Kabul and Indus rivers provide inflows (∼159.3 km 3 yr −1 ) to the IBIS, whereas flows from the Sutlej Ravi and Beas rivers (∼7.8 km 3 yr −1 inflow into Pakistan) are largely diverted to India as per the Indus Water Treaty (UN 1960, IRSA 2020. To mitigate the diversions to India, the Mangla Dam (current capacity of~9.1 km 3 , storing 36% of the total inflows) and Tarbela Dam (current capacity of~7.4 km 3 , storing 12% of the total inflows), 16 barrages, 14 interriver link canals and ∼60 000 km of main canals were built to provide controlled irrigation to ∼16 million ha (FAO 2012, Stewart et al 2018. For a detailed schematic diagram of the IBIS (depicting the system as of 2020) refer to figure A1 in Ahmad et al (2020a).
Surface water resources in the IBIS are allocated seasonally at the provincial level-combining a complex procedure of water resource assessment and forecasting-and then to canal command areas grouped in management zones (MZs) (GOP 1991). The irrigation system is supply-based, originally designed to support one crop per year while achieving equity of surface water supply between users proportional to their landholding (Ahmad et al 2009). With groundwater development over the last five decades (mainly in the Punjab province) (Kirby et al 2017), irrigated area and crop intensity increased (Ahmad et al 2014a). Water conveyed to farmers through an extensive canal network follows a time-based (generally seven-day) rotational system, where in order to achieve equity in surface water supplies (Bandaragoda 1998) each farmer can take the entire flow of the outlet for a period proportional to their land area and not by crop demand. Depending on the availability of water resources, this commonly results in seasonal differences in the cropped area and/or deficit irrigation practices (Ahmad et al 2019). In times of water shortage, farmers either apply deficit irrigation and/or use groundwater (where possible).
Assessing irrigation performance in large systems like the IBIS is complex because of the gaps in water supply and demand information (Bhatti et al 2019, Bierkens and Wada 2019). Several performance metrics have been proposed to evaluate irrigation systems (Burt et al 1997, Karimi et al 2012. Equity is the key performance indicator used in the IBIS management. Traditionally, equity is calculated from canal water supply (i.e. cumecs per hectare), but flow monitoring stations and associated data in the IBIS are not generally available for all irrigation outlets, have limited geographic spread, or require laborious compilation and harmonization by various line agencies (Anwar and Bhatti 2018). In a water-scarce system like the IBIS, equity in water use is more relevant for farmers and local to national management and can be computed from satellite-driven actual evapotranspiration (ET a ). Satellite-driven ET a has long been recognised as a robust and cost-efficient approach to quantify overall water use and monitor water management practices (Bastiaanssen and Bos 1999). By including all possible sources of water depletion (including irrigation supply, water from rainfall and soil water) in a spatially-explicit way, satellite-driven ET a provides a more complete picture about the state of water resources in an irrigated area than canal water supply alone (Bastiaanssen and Bos 1999). Satellitedriven ET a was used in a limited number of canal commands in the IBIS for assessing irrigation performance (Ahmad et al 2009, Usman et al 2015, Shah et al 2016. ET a characteristics were assessed in the IBIS (or large parts thereof) using locally adjusted satellite-driven models (Bastiaanssen et al 2002, 2012, Ahmad et al 2008, Liaqat et al 2015, Simons et al 2020. Except for the 8 year data period used in Simons et al (2016) to assess the long-term water balance components in the system, the other cited studies estimated ET a for a particular year only. The IBIS has large spatial and temporal precipitation variability (Cheema and Bastiaanssen 2012), high and variable evaporative demand (Bastiaanssen et al 2002, Ahmad et al 2019, high seasonal variability in inflows (Archer 2003), variability of the naturally saline soils (Qureshi and Barrett-Lennard 1998) and variability in groundwater quality and use (Qureshi et al 2014, Ahmad et al 2014b, which all lead to variable cropping systems (Cheema et al 2020). Besides these dynamic atmospheric and landscape processes, infrastructure projects enhancing the IBIS's irrigation network, regulation and management continue to be implemented (Young et al 2019). To this day, no study has attempted to describe the time-series inter-seasonal ET a or irrigation performance dynamics of the IBIS. The aim of this letter is, for the first time, to assess two metrics of irrigation performance, seasonal irrigation variability and equity in the early 21st century (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) for the 46 IBIS operational canal commands downstream of the Tarbela and Mangla dams (figure 1). Verified satellite-driven ET a estimates are used for this purpose. The assessment is conducted for the two main irrigation seasons: wet summer 'Kharif ' (April to September) and dry winter 'Rabi' (October to March). The approach utilises only the characteristics of the probability density function (PDF) of the ET a pixels within the canal command to assess the irrigation performance metrics.  Two irrigation performance metrics, seasonal variability and equity were assessed using ET a as a proxy for water use based on the concepts of Bastiaanssen et al (1996) and Siddiqi et al (2018): (a) Seasonal variability is defined as the degree of temporal consistency in seasonal inter-and intra-canal command water use, and can be inferred from past seasonal estimates of ET a . (b) Equity is defined as the inter-and intra-canal command variability of water use for a particular season and can be inferred by the spatial variability of the ET a estimates within the boundaries of the canal command. It is noted that ET a does not equate to canal deliveries as water can be lost to recharge and return flows (Ahmad et al 2014a). Also, ET a can be supplemented from groundwater extractions in some canal commands with access to good quality groundwater (Ahmad et al 2005). Using only ET a derived metrics cannot address the adequacy of water supply to satisfy crop water demand (Moran 1994). However, in a system in which equity is the primary management aim and deficit irrigation is commonplace; the proposed performance metrics are suitable to assess irrigation performance and inform irrigation management.

Materials and methods
Evaluating ET a within canal command areas in time and space therefore can provide a way to diagnose seasonal variability and equity. Time-series of violin-plots (Hintze and Nelson 1998) are used for this purpose, alongside the coefficient of variation (CV, the ratio of the standard deviation to the mean expressed here as a percentage). Violin-plots combine traditional summary statistics of box plots and the shape of the PDF in a single plot, making them useful to assess the dynamics of PDFs both in space and time. The violin-plots are used to assess the intra-and interseasonal variability of a canal command, all pixels within a canal command boundary are included to capture the inherent variability of the existing land cover in space and time. The variability in a given season may include expanding or contracting irrigated fields (depending on the water resource availability or construction/abandonment of canal infrastructure) or encroachment of semi or urban areas on agricultural land. The assessment is largely independent of the accuracy of the ET a estimates if these are sensitive to the seasonal variations and magnitudes of the different land covers present in the IBIS, and do not show systematic biases. The verification (see supplementary material) of the IBIS CMRSET ET a estimates (Peña-Arancibia et al 2020a) used herein confirmed they are fit-for-purpose.
The seasonal non-parametric Mann-Kendall test (Hirsch et al 1982) and the Sen-slope method (Sen 1968) were used to investigate trends in the seasonal ET a time-series, both per-pixel and aggregated to the canal command scale.

Diagnosing seasonal variability and equity using violin-plots
Three representative canal commands areas, in locations with different irrigation management, were chosen to illustrate how seasonal variability and equity can be interpreted using violin-plots from the ET a PDFs (figure 2). The three canal commands have long distances (>20 km) between head-end and tailend irrigators to test equity. The season chosen is the wet Kharif, when high water demand crops (e.g. rice and cotton) are grown. The first canal command assessed, the Upper Chenab Canal-Internal (UCC-INT) canal command (area 3 in figure 1), is in the upper IBIS, precipitation is ∼520 mm yr −1 , surface water irrigation deficits are normally met by pumping groundwater (Ahmad et al 2009(Ahmad et al , 2014b. Rice and wheat are the major crops in the Kharif and Rabi seasons, respectively. The PDF of Kharif ET a shows that generally ET a > 400 mm season −1 , with most pixels having an ET a > 500 mm season −1 . The intra-canal command CV is low at 8.95%. The violin plot shows that UCC-INT can be classified as a reasonably equitable canal command, if both surface and groundwater are considered. Variability can be explained by the temporal changes reported in the time-series of Kharif ET a violin-plots. If the shape and magnitudes are relatively maintained through time, then UCC-INT can be classified as having low seasonal variability. The second canal command assessed, the Abbasia canal command (area 30 in figure 1), is in the middle IBIS, precipitation is ∼189 mm yr −1 , there is little groundwater pumping to supplement surface water supplies due to marginal to poor groundwater quality (Ahmad et al 2014b). Cotton and wheat are the major crops in the Kharif and Rabi seasons, respectively. The PDF of Kharif ET a shows reasonably high ET a variability. Very low ET a (<100 mm season −1 ) can be attributed to the desert areas to the south of the canal command. However, tail-end areas covered by the irrigation network also show low (<300 mm season −1 ) ET a . Most pixels have ET a > 400 mm season −1 . The intra-canal command CV is high at 49.59%. Although some of the CV can be explained by the desert areas, the Abbasia canal command can be classified as reasonably inequitable.
The third canal command assessed, the Pinyari canal command (area 46 in figure 1), is at the IBIS system-end, precipitation is ∼154 mm yr −1 , groundwater is of poor quality and therefore there is little or no pumping to supplement surface water supplies (Ahmad et al 2014b). Rice and wheat are the major crops in the Kharif and Rabi seasons, respectively. The PDF of Kharif ET a shows high ET a variability. Very low ET a (<100 mm season −1 ) can be attributed to the desert areas south of the canal command. However, many of the tail-end areas covered by the irrigation network show low (<300 mm season −1 ) ET a . ET a is highly variable, and in most pixels, it generally ranges from 200-450 mm season −1 , low when compared to the other two canal commands. The intra-canal command CV is high at 45.18%. Note that the Pinyari canal command CV is lower than the Abbasia canal command, but it is evident that because of the variability shown in the shape of the violin plot (i.e. the large widths with ET a in the range of 300-400 mm season −1 in figure 2(i) compared with those for this ET a range in figure 2(f)) it is overall more inequitable. This occurs because the distributional characteristics between canal commands vary and CV is not a robust estimator when data are not distributed normally (Ospina and Marmolejo-Ramos 2019). Figure 3 shows the spatial characteristics of mean annual (cropping year, April to next March), mean Kharif and mean Rabi ET a , respectively. Figure 3 also shows spatial ET a trends (computed using Senslope per pixel) and ET a CV. Spatial differences in mean ET a are more noticeable during the wet Kharif season, with relatively higher values in the Punjab mixed zone serviced by the Sidhnai Barrage (∼500 mm season −1 , orange hues in figure 1, see table S1 for details), and in the Sindh rice area serviced by the Guddu Barrage (∼530 mm season −1 , dark blue hues in figure 1). The lower values are at the system-end canal commands serviced by the Kotri Barrage (∼395 mm season −1 , green hues in figure 1). Dry Rabi ET a values are less variable, with relatively higher values in the Sindh rice zone serviced by the Sukkur Barrage (∼430 mm season −1 , red hues in

IBIS ET a spatial and temporal variability and equity
the probability density function (PDF) of the canal command ETa pixels and a non-parametric (kernel) fitting and (c) the violin-plot using the non-parametric (kernel) fitting and associated box and whisker plot; (d)-(f) idem but for the Abbasia canal command; (g)-(i) idem but for the Pinyari canal command. figure 1). Trends in both seasons seem to be generally slightly increasing (5-10 mm season −1 yr −1 ), except for the high increases (>10 mm season −1 yr −1 ) in the Chashma Right Bank Canal in Khyber Pakhtunkhwa and Punjab provinces during Kharif (areas 18 and 19 in figure 1). There are also decreases (5-10 mm season −1 yr −1 ) in the Sindh rice areas served by the Guddu and Sukkur barrages. The CV shows mostly low values (<20%) in areas covered by the irrigation canal network.
The Kharif ET a violin-plots time-series for the 46 canal commands are shown in figure 4. The statistics reported for each time-series are the mean annual CV (in %), the Sen-slope trend (in mm season −1 yr −1 , computed as the trend of the spatial mean ET a values) and the trend p value. The canal commands in the third and fourth row in green correspond to the Punjab SVC MZ and are used as an illustrative example because they show some variation in the PDF shapes and some variability. ET a is increasing significantly (p < 0.01) in all canal commands (mean of 4.8 mm season −1 yr −1 range 3.37-8.96 mm season −1 yr −1 ). CVs are low on average (mean of 15.70%, range 9.70%-33.20%). The timeseries mean ET a is noticeably increasing in Karanga, LDC, UPC + UMC, FC and UBC + QC (areas 10, 12, 13 and 14 in figure 1). These canal commands can be considered with low variability and equitable. The PDF in ESC (area 15 in figure 1) shows more variability on a season-to-season basis, this canal command is at the edge of the system and adjacent to desert areas with irrigation expanding therein (Kirby and Ahmad 2016), as the PDF shape in later seasons has generally higher ET a values in the lower range (see the areas with low ET a in figure 3(d) alongside the increasing ET a trends in figure 3(e) for this canal command). Other canal commands in different MZs also share ESC's characteristics (the Thal canal command, LBC, Abbasia, PFC, Ghotki, NWC, Dadu, U Nara, respectively areas 17, 29, 30, 32, 36, 37, 39, 43 in figure 1). These canal commands can be considered with reasonable variability and reasonably inequitable. CRBC-P (area 19 in figure 1) is a special case, it shows a steady increase in mean and PDF ET a , becoming more stable in later seasons. CRBC-P was completed in 2003 thus expanding the irrigation frontier (Atta-Ur-Rahman and Khan 2008).
It is readily apparent that in the Punjab FLC MZ (two first rows in blue hues), the canal commands with access to groundwater to supplement surface water shortfalls have less variability and are relatively equitable. The PDFs are generally narrowly concentrated around the median, and median ET a values do not vary much seasonally, CVs are low on average (mean of 10.90%, range 7.60%-15.16%). Some canal commands particularly in the Indus mixed MZ (orange hues, see particularly URC, Haveli, Sidhnai) also have these characteristics, they have access to groundwater and inter-river transfers from the Indus River (Ahmad et al 2014b). Canal commands in other MZs further downstream, generally show PDFs that are oblong around the median, which means more intra-canal command variability in ET a and hence relatively less equity. This can be seen in all the canal commands in the Sindh Guddu Barrage MZ (dark blue hues), Sindh Sukkur Barrage MZ (red hues). There are no drastic changes in the shapes of the canal command PDFs or drastic  changes in the median values; hence although relatively less equitable, the inter-season variation is low. In the Sindh Kotri and Kalri barrages MZ, the Kalri and Lined canal commands show very similar PDFs and medians over the years. Although ET a means are relatively lower, these canals can be considered of low variability and relatively equitable. Conversely, Pinyari and Fuleli show changes in both PDFs and ET a means over the years, thus these canal commands can be considered of high variability and relatively less equitable. There are significant (decreasing) increasing trends (p < 0.01) in ET a in (2) 32 canal commands.
The Rabi violin-plots time-series for the 46 canal commands (figure 5) show similar patterns as the Kharif violin-plots, but overall, there seems to be more intra-canal command equity as PDFs do not change drastically (even for Pinyari and Fuleli). There are less variations, and slight significant increasing trends (p < 0.01) in ET a in 24 canal commands.

Discussion
This study used satellite-driven ET a to provide a consistent spatiotemporal way to assess two commonly used irrigation performance metrics, seasonal variability and equity, in the 46 canal commands in the IBIS. The main advantage is arguably that satellite and meteorological data needed to estimate ET a are readily available online or through the developers, or more recently via Google Earth Engine (e.g. Zhang et al 2019). Monthly CMRSET ET a data from March 2000 to December 2018 for the IBIS area subtending the canal commands (figure 1) at 500 m resolution are available to the public via the CSIRO Data Access Portal (see data availability statement) https://doi.org/10.25919/18h2-vf51. The ET a data was used to estimate the time-series water balance for the canal commands in the IBIS , thus providing information on magnitudes and trends of rainfall, ET a , canal deliveries, groundwater levels and unaccounted inflows and outflows.
Satellite-driven methods that rely solely on routinely used vegetation indices to estimate ET a have known limitations (Glenn et al 2010). Satellitedriven ET a estimates are most accurate for homogeneous land covers. At 500 m resolution, a MODIS pixel is much larger than the average irrigation field in the IBIS, thus the large land cover heterogeneity is expected to result in less ET a accuracy (Glenn et al 2011). The time-lag between the onset of vegetation moisture stress and its manifestation through vegetation indices can result in higher ET estimates than those measured, for example, by complex soilvegetation-atmosphere transfer models (Kustas and Anderson 2009). In addition, the accuracy of satellitedriven ET a models is related to the accuracy of the ground methods on which they are calibrated (Glenn et al 2011). CMRSET was calibrated with flux tower data, which have an accuracy of 10%-30%. Under semi-arid conditions in Australia (similar to those in the IBIS), CMRSET ET a estimates had comparable or better accuracy than other satellite-driven ET a models, when compared to flux tower data (Guerschman et al 2009). An independent evaluation of CMRSET using long-term catchment water balance data from 227 Australian catchments also showed an accuracy (in terms of the root mean squared error) of ∼20% with no significant biases (Guerschman et al 2009). Regardless of the uncertainty in CMRSET's ET a accuracy, the nature of the present assessment relies on the inter-seasonal patterns of ET a within and between canal commands. If these ET a estimates are sensitive to the seasonal variations and magnitudes of the different land covers present in the IBIS (as described in the supplementary information), they provide reliable understanding of canal command seasonal variability and equity, and hence are useful for irrigation management. Although the 500 m MODIS spatial resolution may be insufficient to assess irrigation performance for individual irrigation fields, it is suitable at the canal command scale. It is possible to improve the spatial resolution for areas that require more detail using data from the recently launched satellites like Landsat 8 (30 m) and the and Sentinel-2 (10 m) satellite constellations, alongside blending techniques that make use of mediumresolution-high-frequency MODIS data with highresolution-low-frequency Landsat/Sentinel 2 data (Emelyanova et al 2012). The method can also be used to assess if targeted (non)structural irrigation management has any influence on ET a variability and equity. This assessment of the management of the nearly fully utilized IBIS (Kirby et al 2017, FAO 2019 revealed that it experienced low inter-seasonal variability but with inequity (to varying degrees) in canal commands that did not have ready access to groundwater (particularly in Sindh). High variability and relatively low equity were found in canal commands at the system-end, mainly during Kharif. This result concurs with Hassan et al (2019), that analysed reliability in Sindh using flows in the Sindh barrages. It was found that ET a is increasing in most canal commands during Kharif. Temperature, the basis of the Hargreaves reference ET 0 used herein in the absence of more spatially resolved data, may not be the strongest control of ET a as indicated by experimental studies elsewhere (Donohue et al 2010, Hobbins et al 2012. Moreover, although 19 years of ET a estimates were used in the assessment, thereby capturing this century's management rules, climate variability is captured only to an extent (Cook et al 2013).
Irrigated agriculture has expanded in recent decades in Punjab, chiefly because of groundwater pumping (Ahmad et al 2014b(Ahmad et al , 2019, whilst canal command water deliveries remained essentially constant (Ahmad et al 2020c). Although our study suggest that groundwater use has achieved low variability and equity in most Punjab canal commands, there are concerns about the sustainability of groundwater because of large decreases in groundwater tables in some canal commands (Watto and Mugera 2016). Recent studies (Young et al 2019, Simons et al 2020 suggest groundwater withdrawals are compensated somewhat by recharge, but there is evidence from GRACE total water storage anomalies and hydrological modelling that this balance maybe tilting to an unsustainable situation (Wada et al 2012, Iqbal et al 2017. Water balance modelling for the IBIS by Ahmad et al (2020a), which used long-term observations on groundwater levels, confirmed that water use is becoming less sustainable in several canal commands, especially in the Punjab province. The challenges are further exacerbated by changes in the seasonality and amounts of rainfall, and by rapidly reducing current limited (∼10% of mean annual flows) water storage capacity due to sedimentation (Ahmad et al , 2020c.

Conclusion
This study showcased how long-term satellite-driven ET a , and a visual representation summarizing the inter-and intra-canal command ET a pixel variability through violin-plots, can be used as a pragmatic proxy to infer irrigation variability and equity in large irrigation systems, two metrics commonly used to asses irrigation performance. The assessment conducted per irrigation season (wet 'Kharif ' and dry 'Rabi') in the 46 IBIS canal commands areas (∼160 000 km 2 ) from 2000 to 2018, shows clear patterns: (a) canal commands that supplement canal water supply with groundwater are reasonably of low variability and equitable (mostly in the Punjab province); (b) canal commands in most MZs that rely mostly on canal water supply have low variability over the years, and show intra-canal command inequities; (c) canal commands at the system-end are outstanding in that they can have high variability and also show variable intra-canal command inequities. While good quality groundwater can make up for shortfalls in canal water supply in some canal commands and improve equity, its use over time can become unsustainable and hamper irrigation management goals. The use of satellite-driven ET a time-series provides an opportunity for monitoring irrigation dynamics and the assessment of structural and policy improvements on irrigation performance, in the IBIS and irrigated systems elsewhere.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.25919/18h2-vf51 (Peña Arancibia and Ahmad 2020.