Estimation of glacier mass balance using remote sensing and GIS technology in the Hindu Kush region of northern Pakistan

The Hindu Kush range of northern Pakistan has a large variety of Alpine and valley glaciers. These glaciers are subjected to volumetric changes owing to a host of factors (i.e. altitude, slope, glacier types, and aspect). The geodetic survey and field-based data are not available for assessment of glacier’s Equilibrium Line Altitude (ELA) for Hindu Kush due to rough topography and rigorous climate. However, the mass balance of the alpine and valley glaciers can be obtained through ELA and Accumulation Area Ratio (AAR). Therefore, this study deals with the estimation of ELA, the hypsometrically controlled methods based on AAR and Accumulation Area Balance Ratio (AABR). Various ratios between 0.4–0.8 with 0.05 interval for AAR and 0.67–4.4 with 0.01 for AABR are applied to obtain the ELA of valley glaciers. Moreover, the glacier contours with the interval of 50 m ± 5 m were employed for these estimations. The constant AAR and AABR ratios rather than constant glacier area adjust glacier hypsometry based on ELA variations. The results show that for AAR (0.5), 50 m ELA declined from 5409 to 5359 m and reduced (4%) the geometric area. The maximum decrease of 150 m ELA for these glaciers is reported for the AAR-AABR ratio of 0.1–4.4, showing a significant decrease. The estimated decrease in the glacier area of Hindu Kush is extremely variable specifically in Phargam, Khorabohr and Gordoghan.


Introduction
Northern Pakistan is covered with the world's highest mountain ranges including the Karakorum, Himalaya and Hindu Kush. These ranges stretch over 128,731 km 2 in northern Pakistan. About 28% of the Hindu Kush, 71% of the Karakoram, and 34% of the Himalayan areas lie above 4000 m elevation (Ashraf et al. 2017). Elevation of these ranges varies between 1000 m in the south and 8500 m toward the north and northeast. At high elevations (up to 300 m and above) these ranges are covered with glaciers. Altogether 5218 glaciers have already been identified which cover a total glaciated area of 15040 km 2 . These glaciers are considered water /tanks for the sub-continent, providing water for various uses (Cogley 2011). These glaciers contribute a total ice reserve of about 2738 km 3 (Ashraf et al. 2017;Gilani 2015).
Snow/ice and glaciers are most sensitive to climate change and global warming. It is reported that the glacier's extent is shrinking faster than ever, especially at the lower elevation in the northern Pakistan mountainous ranges (Ashraf et al. 2017;Rasul et al. 2008), which may create water scarcity in coming years. Over the recent years, the glaciers on the Tibetan Plateau and its surrounding mountain in the Himalayas and the Hindu Kush have been retreating faster, whereas the glacier in Karakoram Range are either stabilized or re-advanced in nature (Rasul et al. 2008).
Globally, an increasing number of glaciological studies are focusing on glacier monitoring and glacier changes extent, mass balancing in areas such as Alaska, Patagonia, Andes, Alps, Himalaya and Central Asia. Majority of these studies used various techniques based on satellite images and Digital Elevation Model (DEM) (Racoviteanu et al. 2008;Brun et al. 2017;Konig et al. 2001;Nussbaumer et al. 2001;Trondheim 2013;Zemp et al. 2013). Although, the field-based survey for the estimation and assessment of ELA (ELA, contour line separating the area between accumulation and ablation) is not available for the Hindu Kush range due to harsh climate and rugged terrain especially above 5000 m (Rashid et al. 2015). Therefore, application of indirect methods based on remote sensing techniques are widely applied to determine the ELA and AAR for the HKH glaciers. The AAR-AABR is an extensively useful technique for the estimation of ELA and net mass balance over one year (Benn and Ballantyne 2005;Benn and Lehmkuhl 2000;Lukas 2007;Pellitero et al. 2015). The ELA is the zone with elevation and where the accumulation and ablation are equal over one year, so the mass balance compared to elevation is zero (Cogley 2011). The surface area covered with the glacier away from the ELA (either positive or negative) and has a greater influence over the glacier mass balance as compared to the area near to ELA (Baig and Muneeb 2021). For the constant glacier area the ELA is equal to the zero net balance, whereas it is highly varied for the annual changes in the ablation and accumulation (Das and Rai 2018;Pellitero et al. 2015). Usually, ELAs corresponded with the local topography and climate change. Thus, the changes in the climate past and present-day can be monitored and understood from ELA (Donghui et al. 2017;Wu et al. 2020).
In the pretext of the above literature, this study was designed to estimate the ELA-AAR and AABR and hypsometry for the glaciers of Buni Zom valley of the Chitral River Basin (CRB), in the Hindu Kush region of northern Pakistan. The AAR-ELA techniques were employed for the estimation of glacier mass balance from multispectral data and digital elevation models.

Study area
The CRB hosts 806 large and small size glaciers (Shakir et al. 2010),covering an aggregated area of 1815.844 km 2 , which accounts for 12.3% area of the whole basin (SUPARCO 2017). This study focuses on three glaciers (Gordoghan, Khorabohr and Phargam) lying on the eastern edge of the river basin, having an area greater than 15 km 2 (Fig. 1). These glaciers are part of Booni Zom group range, which is one of the famously known mountain range having the highest peak namely Booni Zom (main) with an elevation 6542 m (21463 ft) (Harlin 2008). The selected glaciers are situated between 72° 15′ 29'' to 72° 23′ 29'' E Longitude, and 36° 07′16'' to 36°11′48'' N Latitude with an elevation from 3700 to 6200 m Above Sea Level (ASL) on the eastern side of the CRB (Fig. 2). Along with other glaciers, the selected glaciers feeds the Chitral River from spring to summer season of the years. Originating from the westerlies, the maximum annual precipitation (730.3 mm) is recorded in the form of snowfalls in the winters and spring seasons . Moreover, the mean monthly temperature of the coldest months (December to January) remain 8 ºC and the hottest month (June to August) is 28 ºC is recorded. Glacier melt from the watershed locality contributes significantly to River Chitral as it runs down the valley (Ahmed 2020). Tirich Mir, lie in the catchment area of the CRB, with a mean elevation of 7708 m and is the highest snow peak Hindu Kush region (Gul et al. 2020).

Glacier data
The temporal satellite images (Landsat-5TM, 7ETM + & 8OLI) were acquired from the U.S. Geological Survey (http:// earth explo rer. usgs. gov/) to map the glacier's outline of selected glaciers (Gordoghan, Khorabohr, and Phargam glacier) for the year 2018 using standard band combination and Normalized Difference Snow Index (NDSI) methods in Arc-GIS (version 10.5). From the ASF website the Advanced Land Observing Satellite (ALOS) Phased Array Type L-band SAR (PALSAR) 12.5 m resolution DEM was downloaded from the https:// earth data. nasa. gov/ eosdis/ daacs/ asf, to calculate the hypsometry of the study area and to extract the contours for the ELA and AAR estimation. The updated glaciers inventory of the whole CRB and Buni Zom valley was downloaded from Global Land Ice Measurements from Space (GLIMS) available on www. color ado. edu (GLIMS and NSIDC 2018).
The AAR and ELA were generated using digital image processing techniques to estimate the glaciers mass balance of the three selected glaciers of the Buni Zom valley.

Zero mass-balance estimation of AAR and AABR
In the present study, the hypsometrically controlled AAR and AABR methods were employed for ELA estimation of mass balance and glacier climate reconstructions were applied in the Buni Zom area (CRB) of the Hindu Kush region.
The over-glacier area is divided into two zones by the equilibrium line, the accumulation zone and the ablation zone. The accumulation area (ice deposit) is measured by the end of the hydrological year, separated by the equilibrium line from the ablation area. In this study, AAR is the ratio between the accumulation area divided by total area of the entire glacier (Eq. 1). In the ablation season i.e. September, the glacier maximally exposed to climate changes (temperature). The characteristic of the glacier surface in the accumulation zone is separated from the ablation zone by using satellite imagery allowing possibility to map the ELA.
Based on the Kurowski method (Eq. 1), a sophisticated toolbox for processing of glacier ELA based on AABR and AAR was proposed by Pellitero et al. 2015 in the ArcGIS platform. To run the toolbox in ArcGIS the Python code has been used on the glacier surface reconstructed from DEM as an input. The process of ELA determination becomes (1) AAR = Accumulation Area∕Total Area of Glacier easy and simplified through this toolbox. This toolbox can be effectively applied over a single or large dataset of several glaciers (Pellitero et al. 2015). In this study, the tool proposed by Pellitero et al. 2015, was used to compute the ELA-AAR and ELA-AABR of the selected glaciers.

Estimation of ELA-AAR and ELA-AABR for glacier mass balance
Various size and valley type glaciers are found in the Hindu Kush region. The mass balance of the valley or Alpine-type glaciers was determined by the ELA (Kern and Laszlo 2010;Pellitero et al. 2015).
The mass balance of the glaciers was measured through the ELA and AAR techniques. The glacier area was divided into the accumulation and ablation zone by the ELA line. The accumulation zone lies in the upper part of the ELA line and ablation zone lie below the ELA line of the glacier, where the loss of mass is noticeable. Glacier mass balance is acquired with ELA techniques to assess the hypsometric control of glaciers. An automatic tool developed by Pellitory  , (2015) was used to calculate the ELA of these glaciers (Gordoghan, Khorabohr, and Phargam) of Buni Zom from the given AAR and AABR value. This procedure satisfied the glacier mass balance for different altitude glaciers using DEM. In this study, the lowest standard deviation value for the AAR and AABR is 0.5 and 0.6 were employed to extract the ELA, as these three glaciers of CRB lie slightly in different climate zone (Baig and Muneeb 2021). Initially, zero net balance ELAs was assumed for each glacier of the CRB and assumed these are in equilibrium with the climatic condition. Globally, for the AAR and AABR with the lowest standard deviation, the eleven evenly spaced ratios from 0-1 (0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0) were applied correspondingly with altitude (Kern and Laszlo 2010;Pellitero et al. 2015;Baig and Muneeb 2021). Depending on the location, size and types of the glaciers the studies of Rashid (2015), Igneczi and Nagy (2013), Kern and Laszlo (2010) have proposed various AAR and AABR ration for the estimation of the ELA (Table 1).
For the study glaciers AAR-AABR for several ratios ("0.4 0.8 0.05" "0.9 4.4 0.1") were calculated and AAR-AABR was 5409 and 5359 m. For each glacier the DEM were used for the balance contour intervals and the resulting ratio of ELA-AAR and AABR were stored (Fig. 3). Based on the above technique the mass balance of the glaciers was used to calculate with the help of AAR-AABR. The glacier mass balance has been influenced greatly (either positively or negatively) by the altitudinal area away from the ELA as compared to an area that exists adjacent to the ELA. For the estimation of the AABR the toolbox follows Osmaston's (2005) technique and for the AAR calculation, it follows the González- Trueba and Serrano (2004) approach (Osmaston 2005;Trueba and Canadas 2004). The outputs show a constant ratio for the accumulation and the ablation area because the glaciers of the Buni Zom are in a steady state (0.5).

Hypsometry of the glaciers
As given in Fig. 4 the Gordoghan the largest glacier occupying 17.22 km 2 ± 0.21, followed by Khorabohr at 14.98 km 2 ± 0.25 and Phargam at 9.88 km 2 ± 0.84. These glaciers have a continuous-time series change in the area during the ablation season and accumulation season. The studied glaciers lie in the extreme altitudinal range of the valley. The largest glacier Gordoghan of the study lies between 3714 and 6193 m with an area of 17.22 km 2 , followed by the Khorabohr glacier lies between 4328 and 6063 m with an area of 15 km 2 and the smallest glacier Phargam between 4225 and 6110 m with an area of 11 km 2 . These glaciers were mostly clean and ice-covered. Whereas, the highest mean elevation of 5177 m (ASL) was observed for the Phargam glacier, followed by the Khorabohr glacier (5195 m (ASL)), and Gordoghan glacier (4953 m ASL) as given in Fig. 4. These glaciers significantly contribute to melt water in the Mastuj and Chitral rivers and indicated the continuous shrinking of the glaciers. The detail of more than 250 glaciers of different sizes, locations, and elevations of the CRB including Buni zoom valley is given in Fig. 1

Zero glacier mass balance
The zero net balance was estimated for the ELAs and expected that these three glaciers of Buni Zom valley are in equilibrium with climate conditions. The corresponding net-balance ELAs were estimated from the reconstructed glacier area and contours for the three studied glaciers of Buni Zom valley. The recommended value AAR (0.5) for the valley glaciers, the linear regression was applied over the ELA (5359 m) values in this study (Fig. 5). Based on the AAR-ELA method the negative or positive gain of mass/ mass balance in the Buni Zom valley was found through the linear regression model based on the Correlation coefficient (Eq. 2). The results show a negative net balance (-0.4806) for AAR-ELA for mass balance. Tawde and Kulkarni (2016) found a good correlation between the AAR measure and field mass balance of the Baspa basin glaciers (Shaune Garang and Gor Garang) of Western Himalaya and developed a correlation among these two variables (Tawde and Kulkarni 2016).  (2017), map (b) Buni Zom valley, glaciers with smallest sizes (< 1 km 2 ) are shown with smaller pink dots, while largest sizes between 60 and < 80 km 2 with larger pink dots. Glacier outlines and their hypsometric curves were derived from ALOS-PAL-SAR DEM data

Glacier mass balance based on ELA-AAR and AABR
Due to rough topography and various types of glaciers (small to large size), the ratios between 0.4 and 0.8 with intervals of 0.05 for AAR-ELA and AABR-ELAs between 0.67 and 1.7 with 0.01 intervals were used to calculate the ELA for the studied glaciers. Along with the above parameter, the glacier contours with the interval of 50 m ± 5 m error default were selected to calculate the different ELA-AAR. To adjust the glacier's geometries (area, extent), the constant AAR ratio approach was used to calculate the ELA. The constant AAR of 0.5 with 50 m ELA declined from 5409 to 5359 m of selected glaciers to adjust the glacier geometries by reducing 4% area. The highest reduction  (2013) Outlet 0.58 Pellitero et al. (2015) Valdenievas palaeoglacier in N Spain 0.6 Baig et al. (2021) Karakorum range glaciers 0.6 Kern and Lazlo (2010) 0.1-1 km2 extension 0.44 ± 0.07 1-4 km2 extension 0.54 ± 0.07 Larger than 4 km2 0.64 ± 0.04 All glaciers (World Glacier Inventory) 0.559 ± 0.09

Fig. 3 Map of selected studied glaciers and elevation corresponding to glaciers ELA-AAR & AABR. A minor variation was estimated in AAR and AABR
1 3 The hypsometric curve of selected glaciers (35.5%) was much flatter because they lie at transitional altitude (4000-6000 m). The accumulation area (25%) lies to the steeper area above at 5359 m and the 75% ablation area lies below when the AAR (0.5) was compared at 5359 m ELA. While for AABR (1.75) ratio 35% accumulation area and 65% ablation area lie of overall selected glacial area. An increase in the AAR from lowest (0.4) to highest (0.8) reduces the 200 m ELA of the studied glacier from 5409 to 5209 m, 30% decrease in the accumulation area. Similarly, the rise in the AABR ratio from lowest (0.9) to highest (4.4) declines the 50 m ELA for the studied glaciers from 5159 to 5109 m and 25% loss in accumulation area in Table 2

Climatic effects on glaciers mass balance
The impact of climate change can be certainly monitored and asses over high altitudes as compared to low altitude mountains (Yang et al. 2011). It is evident from the research and studies that the Alpine and valley glaciers are severely influenced by environmental changes. As the high-altitude areas seem to be warming more, and perhaps faster than the rest of the globe (Donghui et al. 2017;Gilany 2015;Igneczi and Nagy 2013;Wu et al. 2020;Yang et al. 2011). One of the clearest and potentially severe impacts is likely to be the retreat of glaciers, loss of mass in some cases, vanishing of many alpine glaciers throughout the world (Yang et al. 2011), and ultimately contributes to sea-level rise and increase in the frequency of GLOFs, glaciers outburst (Wu et al. 2020).
The mean temperature and total precipitation data over 48 years  were obtained from the regional metrological office in Peshawar and have been analysed to find the linear trend of the local climate condition in the study area as the local climate condition is generally associated with glacier dynamics. The short and long-term changes in climate patterns were explored to understand the present and future climatic trend and their association with detected glacier mass changes. Most studies investigate to identify the impact of global warming on the glaciers of the Himalaya, Karakorum, and Hindukush (HKH) ranges (Rasul et al. 2008;Roohi et al. 2002).
For the changing trend of climate data, temperature and precipitation, the nonparametric test Mann-Kendall trend test with Sen's Slope Estimator was used to estimate the magnitude of the trend for time series data Gul et al. 2020;Hussain et al. 2015) for both variables. The temperature has a slight tendency towards an increase while the annual precipitation shows a decrease. For the annual mean temperature, the value of Mann Kendall trend of 0.285  and Sen's slope trend of 0.038 ºC were found. While for the annual precipitation the resulting values for the Mann Kendall trend 0.006, -0.242 mm Sen's slope trend was observed ( Figs. 6 and 7). Trends of annual precipitation and mean temperature for the CRB were observed as positive for the temperature while negative for the precipitation during the study period. Ahmad et al (2018) and Hussain et al (2015) has studied the seasonal trend of temperature and precipitation for the HKH range and observed a positive increase in mean temperature. Which resultant in the gradual change in snout area, mass balance and shifting of ELA from low to high elevation. In the current study, the same trend of negative gain of mass and a gradual shift of ELA was observed.

Discussion
The current study was based on the glacier's hypsometry and mass balance extracted from ELA-AAR and AABR techniques. The techniques were applied over the three selected glaciers of the Buni Zom Valley of North Chitral of the Hindu Kush range using Landsat images and DEM. Glacier mass balance is contributing and controlling regional hydrology (Bakke and Nesje 2011;McGrath et al. 2017) and sea-level rise as it plays a key role to evaluate the advance and decline of water storage in glaciers (Gardner et al. 2013;Tawde and Kulkarni 2016). The snow and glaciers cover play an important role in the control of regional and local changing climate dynamics (Das and Rai 2018;Rashid et al. 2015), particularly throughout the past century.
AAR and ELA techniques have been applied to estimate the glacier mass balance by several studies in the past and present for the individual and basin glaciers (Trueba and Canadas 2004;Zasadni et al. 2018). Although field surveys and assessments for every glacier based on the ELA-AAR techniques are not possible, as the snowline and ELA for both season ablation and accumulation are generally estimated through satellite images (Tawde and Kulkarni 2016). The ablation and accumulation are equal throughout the year at ELA point and the zero mass balance is observed at this point (Cogley 2011). The main benefit of the ELA-AAR technique is that it can be used to estimate the mass balance of individual glaciers as well as basins glaciers. The exact location of the ELA on satellite images is thoughtprovoking because of cloud cover, temporal gaps, and fresh snowfall. Hence, the highest observed ELA has been usually used for the estimation of AAR, which result in the gain or loss of the glacier mass (Cogley 2011;Gardner et al. 2013).
The studies of Pellitero suggest a 0.58 AAR value for the North Cascade glaciers and Alpine glaciers 0.7 AAR value with zero mass balance (Pellitero et al. 2015). In this study the value AAR 0.5 and value AABR 0.67 were used for the mid-latitude glaciers (Buni zom. The AAR values (0.4-0.8) and AABR values (0.90-4.4) were applied for the ELA estimation and a decreasing trend was observed for these glaciers. The AABR approach is best suited to obtain the ELA for the clean ice-covered glaciers as compared to the AAR (Pellitero et al. 2015). However, the Osmaston studies suggest the value of THAR (Toe-Headwall Area Ratio) and AAR 1.5 and 2.5, respectively for the East African mountain glaciers with different hypsometry (Osmaston 2005).
Various tools and techniques based on satellite images and digital elevation have been used to monitor the changes in glaciers elevation, retreat and mass etc. (Das and Rai 2018;Gardelle et al. 2013). Donghui et al. (2017) used topographic maps, satellite images, DEM and GPS data to monitor the changes in glaciers area, extent and elevation in selected Yanglonghe glacier no. 1 (5Y432A1) and 5 (5Y432A5) of central China during the time span of 1956 to 2007. The area of glaciers 5Y432A1 and 5Y432A5 has decreased by 4.1% and 15.9%, respectively. The average rate of tinning of glaciers were 20.2 ± 11 m (0.4 ± 0.22 maˉ1) and 16.9 ± 11 m (0.33 ± 0.22 maˉ1) in the ablation zone of glaciers 5Y432A1 and 5Y432A5 (Donghui et al. 2017). Moreover, a similar study was conducted by Tawde and Kulkarni in 2016 to estimate ELA by combining satellite images with Situ metrological data and applying the snowmelt model. He applied this method over twelve selected glaciers in the Chandra basin, Western Himalaya. They found a decrease in the mass of glaciers and compared the traditional AAR technique and the modelled estimates for mass balance which are in good agreement with the geodetic method. The modelled collective mass balance is −1.67 ± 0.72 Gt and −0.79 ± 0.34 m w.e.a from 1999 to 2000to 2009, respectively (Tawde and Kulkarni 2016. The ELA for the Kolohai glacier of Lidder Valley in Kashmir Himalaya lies within 4594 m altitude, while AAR is 0.59. The values of ELA and AAR indicate a negative mass balance in the Kolohai glacier. A slight shift of the ELA was observed toward a higher altitude. This was 4205 m with an AAR value of 0.57 ratios in 2001 and in 2006 it was 4594 m with 0.59 AAR and a negative mass balance was observed for the Kolohai glacier (Rashid et al. 2015). Same AAR value (0.56) was applied in the current study and found the ELA 5359 m with a negative (-0.4806) net balance for the study glaciers was recorded.
In contrast to Himalaya and Hindu Kush Glaciers, the Karakorum glaciers show a positive mass balance. Gardelle et al. documented a slight advance in mass (+ 0.09 ± 0.18 m w.e.a maˉ1) of the Batura glacier, whereas it was 100 ± 0.18 m w.e.a maˉ1 for the Gharesa glacier. This could be the main cause for the different durations of surge events (Gardelle et al. 2013). Lin et al. 2017 documented the mass gain of 0.043 ± 0.078 to 0.363 ± 0.065 m (w.e. yearˉ1) during the years 2011 and 2014 in west Kunlun, eastern Pamir, and the northern part of Karakoram. While the western Karakorum showed a near-stable mass balance of -0.020 ± 0.064 m (w.e. yearˉ1). The Upper Tarim Basin in the Karakorum shows a positive mass balance while a decreasing gradient was noted from northeast to southwest (Lin et al. 2017).
According to Wu et al., the average retreating of glaciers in the Hunza Basin was 0.61 ± 0.07 m from 2000 to 2014. Glaciers with an area of 2878.2 km 2 experienced a mean mass balance of 0.04 ± 0.02 m. Generally, at higher elevation the mass balance shows a significant increase. Glaciers experienced a positive mass balance above 5100 m (ASL) and a negative mass balance below this elevation during the study time period in the Hunza basin. The main reason for this phenomenon is the surging behaviour of the glaciers of the Hunza basin (Wu et al. 2020).

Conclusion
This study finds remote sensing techniques AAR, AABR and ELA were appropriate to estimate the glacier mass balance of Hindu-Kush glaciers. The resolution of ALOS PAL-SAR DEM was enough to estimate the glacier mass balance in the chital river basin of the Hindu Kush range. However, the glaciers boundaries were mapped from published data sets and Landsat images. Quantification and estimation of the glacier mass balance through AAR-AABR-ELA produced quality results. In the automatic tool of ArcGIS to acquire glaciers hypsometry and ELA-AAR and AABR value, reduced processing time, increased the efficiency and quality of the results. The estimated glacier geometries by reducing area are extremely varied; the smaller to larger size glaciers Phargam, Khorabohr, and Gordoghan have reported a noticeable decrease in glacier-ice area. The variation in the ELA between AAR and AABR showed that the ELA-AABR method is appropriate to estimate the ELA for clean and small-size glaciers. We used the AAR value (0.1 to 0.5) and AABR value (0.8 to 4.4) to estimate glaciers zero mass balance of the studied glaciers due to high altitude setting. Glaciers in the CRB experienced a slight negative (-0.4806) mass balance during the study year 2018 with an area reduced by 4.93 km 2 . The changing trend of climate data shows that the temperature has a slight tendency towards an increase while the annual precipitation shows a decrease.
Authors contribution Sidra Bibi conceived the presented idea, and collected the data. Sidra Bibi drafted the manuscript with support from Muhammad Shafique and Neelum Ali. Sidra bibi and Neelum Ali contributed to the interpretation and analysis of the data. Liaqat Ali and Rehman Gul provided critical feedback and helped shape the manuscript. Moreover, Shahla Nazneen and Siddique Ullah Baig prepared the graphs and contributed to addressing the issues raised by reviewers and revising of the manuscript.

Data availability
The co-author should be contacted if someone wants to request the data.