Hydro-climatic Characteristics of Yarlung Zangbo River Basin since the Last Glacial Maximum

Shuang LIU, Kaiheng HU*, Weiming LIU, and Paul A. CARLING

1Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China

2Geography & Environment, University of Southampton, Highfield, Southampton, SO17 1BJ, UK

ABSTRACT Global climate changes significantly impact the water condition of big rivers in glacierized high mountains. However,there is a lack of studies on hydrological changes within river basins caused by climate changes over a geological timescale due to the impossibility of direct observations. In this study, we examine the hydro-climatic variation of the Yarlung Zangbo River Basin in the Tibet Plateau since the Last Glacial Maximum (LGM) by combining δ18O proxy records in Indian and Omani caves with the simulated Indian summer monsoon, surface temperature, precipitation, evapotranspiration and runoff via the Community Climate System Model and the reconstructed glacier coverage via the Parallel Ice Sheet Model. The mean river runoff was kept at a low level of 145 billion cubic meters per year until an abrupt increase at a rate of 8.7 million cubic meters per year in the B?lling-Aller?d interval (BA). The annual runoff reached a maximum of 250 billion cubic meters in the early Holocene and then reduced to the current value of 180 billion cubic meters at a rate of 6.4 million cubic meters per year. The low runoff in the LGM and Heinrich Stadial 1 (HS1) is likely attributed to such a small contribution of precipitation to runoff and the large glacier cover. The percentage of precipitation to runoff was only 20%during the LGM and HS1. Comparison of glacier area among different periods indicates that the fastest deglaciation occurred during the late HS1, when nearly 60% of glacier area disappeared in the middle reach, 50% in the upper reach,and 30% in the lower reach. The rapid deglaciation and increasing runoff between the late HS1 and BA may have accelerated widespread ice-dam breaches and led to extreme outburst flood events. Combining local geological proxy records and regional simulations could be a useful approach for the study of paleo-hydrologic variations in big river basins.

Key words: Yarlung Zangbo River Basin, Indian summer monsoon, glacier change, runoff variability, since the LGM

1. Introduction

Hydro-climatic variations drive the evolution of mountain landscapes and significantly influence the environment of human beings (Xuan et al., 2021). Quantitative understanding of prehistorical hydro-climatic change in a river basin is critical for decision-makers to improve water resource management and mitigate extreme conditions. Strong cooling and warming climate events, especially during the Last Glacial Maximum (LGM), the Heinrich Stadial 1 (HS1),B?lling-Aller?d interval (BA), and the Younger Dryas period (YD), significantly influenced the runoff of major rivers (Teller, 1990; Obbink et al., 2010). The highest major river on Earth, the Yarlung Zangbo River (Chinese part of the Brahmaputra River), originates in western Tibet Plateau and is largely controlled by the Indian Summer Monsoon and glacier activity. Its hydrological system is highly sensitive to climate change events but the underlying causes are unknown (Goodbred, 2003; Sijinkumar et al., 2016; Pickering et al., 2019). Hence, the Yarlung Zangbo River is a suitable case to examine the quantitative relationships between paleoclimatic change and river runoff.

Recent decades of climatic and hydrological data for the Yarlung Zangbo River Basin (YZRB) have been intensely investigated. These studies have included annual and decadal examination of the glacial ice mass balance(Yao et al., 2010); hydrological alterations (Chen, 2012;Sang et al., 2016; Xie et al., 2019); and variability in the freezing depth of soil (Liu et al., 2020a). Such studies have contributed to understanding hydrological responses to future climate change (Gao et al., 2019; Liu et al., 2019b), and quantifying the total water yield of the Tibet Plateau (Wang et al.,2021). In the last decade, geological data, such as glaciallydammed paleolake records and reconstructed hydrological proxies, reflect the paleoclimatic changes to which the YZRB has been subjected, notably since the LGM (Zhu et al., 2013; Nishimura et al., 2014). However, the long-term detail of hydro-climatic variation, specific to the YZRB,remains obscure until improved methods are developed.

Paleoclimate studies are mainly based on isotopic information (Sidorchuk et al., 2008; Veit et al., 2017). For example, δ13C and δ18O series from black spruce tree-rings were used to reconstruct the 1800-2009 discharge of the lower Churchill River, Canada (Dinis et al., 2019).However, the local proxy-based inversion is difficult to apply to the regional scale (Joussaume and Taylor, 1995;Duffy et al., 2019). Climate system models, which can incorporate environmental proxies, reproduce continuous spatiotemporal information and be applied widely to paleoclimate-related studies (Cane et al., 2006), such as runoff reconstruction (Lehner et al., 2017). To develop a deeper insight into climate change since the LGM, a transient climate simulation (termed TraCE-21ka) has been conducted by Liu et al.(2009). Subsequent studies have used the model outputs to explore interesting phenomena, such as abrupt climate events, deglacial monsoon change, and deglaciation (Partin et al., 2015; Tierney et al., 2016; Yan et al., 2020).

To develop the first impression of hydro-climatic variation of the YZRB since the LGM, this study has analyzed data from multiple sources, including climate proxies to indicate Indian Summer Monsoon (ISM) strength, simulated climatic and hydrological data, and ice-dammed lake outburst records.

2. Materials and methods

2.1. Study area

The Yarlung Zangbo River Basin lies between 27°43′-31°12′N and 81°42′-96°54′E, with an area of 0.24 million km2(Fig. 1). It represents the upper reaches of the Brahmaputra River, originating from the Gyima Yangzo glacier at an elevation of 5590 m on the northern slope of the Himalayas (Zeng et al., 2018). The topography declines from west to east with a mean slope of 2.6‰ (Li et al.,2009). Significant tributaries include the Nyang Qu River,Lhasa River, Nyang River, and the Parlung Zangbo. The Yarlung Tsangpo River with a length of 2057 km flows through the South Tibet Valley, and is called the Brahmaputra river when it enters the piedmont alluvial plain of the eastern Himalayas. During 1998-2014, the annual mean precipitation is ca. 675 mm (Xie et al., 2019).The outlet of the Yarlung Zangbo River is at Pasighat,where the mean annual water discharge is ca. 1860 × 108m3(Sarma, 2004).

Fig. 1. Study area, including the location of river basins, elevation, and river network of the YZRB, an outlet of the Yarlung Tsangpo River in China (known as Pasighat, marked with a green star), and the borders of the upper,middle, and lower reaches.

2.2. Available data

Climate proxies to indicate the ISM strength, simulated ISM wind speed, precipitation, and evapotranspiration were accessed. Runoff and glacier area were determined from a climate model and an ice model, respectively, and ice-dammed lake outburst records were consulted. All the available data are listed in Table 1 and briefly explained below.

Table 1. Data used in the current study.

2.2.1. Simulated climatic and hydrological data

We downloaded the decadal-mean annual averages of surface temperature, summer wind speed, precipitation, evapotranspiration, and runoff from the simulation of the Transient Climate of the Last 21 000 Years (TraCE-21ka), which contains outputs from a fully coupled, non-accelerated atmosphere-ocean-sea ice-land surface Community Climate System Model simulation from 22 000 years BP (before 1950 A.D.) to 1990 CE (Liu et al., 2009). The land surface model uses the same resolution as the atmosphere model (3.75° latitude-longitude), and each grid box includes a hierarchy of land units, soil columns, and plant types. Individual glaciers, lakes, wetlands, urban areas, and vegetated regions can be specified in the land units. Nine of global grids from CCSM3 output cover the YZRB.

2.2.2. Reconstructed paleo-glacier coverage

Based on a 1-km climate-ice sheet simulation by the Parallel Ice Sheet Model, glaciation history over the Himalayan-Tibetan orogen during the past 20 ka has been revealed (Yan et al., 2018, 2020). A hybrid stress balance model is used to compute ice velocity, with the shallow ice and shelf approximations for vertical deformation and longitudinal stretching, respectively (Winkelmann et al., 2011).The model has proved effective for modeling ice sheets and ice fields in both past and present climates (Yan et al.,2016). In this study, the glacier areas from a transient simulation were used.

2.2.3. Paleoclimatic proxies

The δ18O series derived from proxy records in Mawmluh Cave in northeastern India (Dutt et al., 2015), and the integrated Holocene proxy based on four cave δ18O records for the Sahiya cave in Northern India (5.7-0.0 ka BP), Hoti cave in Northern Oman (9.6-6.0 ka BP), Mawmluh cave in Northeastern India (33.8-5.5 ka BP), and Qunf cave in Southern Oman (10.6-2.7 ka BP and 1.3-0.3 ka BP) (Liu et al.,2020b) are used to show the ISM change. Stronger ISM could bring more heat and moisture into the YZRB, contributing to increased runoff volume.

Published data of ice-dammed lake outburst floods along the Yarlung Zangbo River (Liu et al., 2019a) were used. Lacustrine sediments at Gega have been dated to 41-13 ka BP using optically stimulated luminescence(OSL) and14C dating, indicating a glacial dam duration of at least 28,000 years (Liu et al., 2015). Dates for upper lacustrine deposits at Jiedexiu slightly vary, with an OSL age of 13.5 ± 0.6 ka BP (Zhu et al., 2013), and14C ages of 14.46 ± 0.17 ka BP (Zhu et al., 2014) and 15.36 ± 0.17 ka BP (Han et al., 2017). The paleo-lake records can provide an indirect indication of runoff changes during the last deglaciation.

2.3. Method

The area-based weighted average of nine-grid information makes the most of the variability on the coarse grid.Then, outputs of climatic and hydrological variables of the YZRB were extracted from the global results. It is noteworthy that river runoff response to precipitation might be delayed at the seasonal scale due to glacier or snow storage.However, while considering decadal-means, as in this study, the delay effect can be ignored. Thus, runoff over all the grids in the YZRB can be considered total decadal-mean annual discharge at the outlet. CCSM3 is a global, coupled ocean-atmosphere-sea ice-land surface climate model with no flux adjustment (Collins et al., 2006). The present mean annual discharge at the outlet (Sarma, 2004) was used to correct the systemic bias of the simulated runoff to show the possible paleo-runoff level. Meridional wind speeds in summer,representing the ISM strength that influences the YZRB,were extracted from the CCSM3 grids over the Bay of Bengal close to the entrance of vapor transport pathway to the river basin. The simulated glaciation over the Himalayan-Tibetan orogen is in the form of a 1-km grid. Based on the boundaries of upper reaches, middle reaches, lower reaches, and whole basin of the YZRB (Fig. 1), the number of 1-km grids representing glaciers was summed to define the total glacier area. Then, glacier area ratio, defined as the ratio of the simulated paleo-glacier area to the present one,was calculated for analysis.

According to international standards, important climate phases for the past 22-ka period has been divided into the LGM (22-19 ka BP), the HS1 (19-14.7 ka BP), the BA(14.7-12.9 ka BP), the YD (12.9-11.7 ka BP), the Early Holocene (EH, 11.7-8.2 ka BP), the Middle Holocene (MH,8.2-6 ka BP), and the Late Holocene (LH, after 6 ka BP).Trend analysis and boxplots were used to examine the change and dispersion of climatic and hydrological elements. A higher dispersion degree of runoff data shown in boxplots reflects a more unstable hydrological status during a certain period. The partial regression analysis tool provided by the Computational and Information Systems Laboratory at the National Center for Atmospheric Research (http://www.ncl.ucar.edu/) was used to check the relative importance of precipitation and surface temperature for glacier change. Standardized regression coefficients characterize the partial regression coefficients in units of standard deviation. The larger the absolute value of the coefficient, the greater the effect on the dependent variable. When focusing on hydrological changes at the tens-of-millennial scale from the last deglaciation, precipitation, evapotranspiration from vegetation and soil, and total runoff are the three most important hydrological elements in a river basin. Therefore, the changes in contributions of precipitation to evapotranspiration and runoff are also discussed. In addition, a power spectral density tool (Welch, 1967) in MATLAB(https://www.mathworks.com/help) was adopted to explore the runoff cycles for different climate phases.

3. Results and discussion

3.1. Runoff tendency and periodicity

Simulated runoff shows a marked peak coinciding with the BA period of abrupt ice-dammed lake outburst events(Fig. 2). In the Holocene, the runoff increased at first and then decreased. Three long sub-periods can be identified in the YZRB 22-ka simulated annual runoff series: ca. 22-15 ka BP, stable low runoff (mean growth rate of 0.6 ×106m3yr-1); ca. 15-9.5 ka BP, rapid growth in annual runoff (mean growth rate of 8.7 × 106m3yr-1); and ca. 9.5-0 ka BP, exhibiting an evident decline (mean reduction rate of 6.4 × 106m3yr-1). Decadal-mean annual runoff in the rapid growth period (ca. 2130 × 108m3) was nearly 1.5 times that in the initial stable period (ca. 1450 × 108m3).Although there is a downward trend in the last 9.5 ka, the average annual runoff is still high (ca. 2040 × 108m3). The smallest mean annual runoff (1450 × 108m3) occurred at the LGM, but the lowest runoff occurred in HS1 with ca.1270 × 108m3(Fig.3). The highest runoff reached ca. 2500 ×108m3in the EH. The largest dispersion degree occurred in the BA period, reflecting an unstable climate condition. The next largest dispersion degree was observed in the LH. The smallest dispersion degree was in the LGM. From the HS1 to the BA, the mean annual runoff increased by 38%.

The large peak value of the PSD estimates indicates a possible periodic signal in the time series (Fig. 4). The time period from the BA to the YD is an inflection point as the BA warming made the runoff cycle duration much longer than it had been before; short cycling was almost non-existent. The sudden climate cooling in the YD terminated the cycle with shorter cycles occurring later. An extremely long duration of high runoff means sustained water energy storage behind the huge ice dam (Gega and Jiedexiu ice dams),which may destroy the dam's stability by infiltration and basal erosion. On the contrary, the abundant short periodic signal in the Holocene can help in understanding the pattern of river water change, which could be in favor of water management, as it perhaps explains why this was a key period for the development of human civilization in China, India, and surroundings areas (Srivastava et al., 2017a). These cycles show the runoff rhythm of the Yarlung Zangbo River at different temporal scales, indicating that the total runoff responded to sustained climate warming but with periodic variations.

Fig. 2. Simulated precipitation and runoff anomaly for the YZRB over the past 22 ka.

3.2. ISM controlling influences

The general trend of runoff provides a satisfactory match with the ISM indicators (Figs. 2 and 5). During the LGM, simulated runoff is low and corresponds with reduced δ18O levels. Increased runoff in the BA warming period is also related to peaks in ISM strength; however,annual runoff does not show a clear response to the YD cooling compared with the δ18O series. The simulated summer meridional wind variation generally agrees with the ISM proxies (Fig. 5). Abrupt increases in monsoon intensity occurred in the EH and the BA-YD period. Multi-proxy records from three separate Arabian Sea sediment cores suggest the primacy of two abrupt increases in monsoon intensity, one between 13 and 12.5 ka, and the other between 10 and 9.5 ka(Overpeck et al., 1996), which is in accord with simulated ISM strength change. In the early 10 ka, the wind speed increased and reached a peak in the EH, while it decreased by ca. 1.5 m s-1from the late EH to the present. Wind speed also responded to climate warming and cooling during the BA-YD period when an obvious increase and drop of speed occurred. Some wind speeds in the LH are very close to those that occurred in the LGM, indicating a trend of weakened ISM has gradually emerged during the Holocene.The wet-dry phase change is similar to climatic oscillations that occurred during the last 15 000 years in the Indian subcontinent (Sharma et al., 2004).

Fig. 3. Boxplots of corrective decadal-mean annual runoff during different periods. The box extends from the first (Q1) to third (Q3) quartile values of the data, with a line at the median(Q2). The triangle represents the average value. The whiskers extend from the edges of box to show the range of the data. By default, they extend no more than 1.5(Q3-Q1) from the edges of the box. Outliers are plotted as separate dots.

3.3. Precipitation determining general runoff change

We found that the annual precipitation drives the change in total runoff (Fig. 2) as the Pearson correlation between the two factors is always high (Fig. 6). However,the square of the correlation coefficient for the LGM is only 0.77 (Fig. 6b), a much smaller value than other periods(Figs. 6c-h). It is noteworthy that the huge glacier coverage and temperature-dominated meteorological condition for ice melting likely resulted in runoff mechanisms more complex than suggested by the current results. Unknown hydrothermal processes in dynamic glaciers make the relation between precipitation and runoff more uncertain. For example, as described by Smith et al. (2017), fluvial supraglacial catchment could attenuate the magnitude and timing of runoff delivered to lower reaches. Interactions between air and ice could also change the water-heat cycle between atmosphere and land surface to affect runoff. In addition, glacier advance may foster the formation of ice-dammed lakes,which would regulate the total runoff by reducing the fluidity of surface runoff.

In the YZRB, precipitation contributes more to the evapotranspiration than the runoff (Fig.7); antiphase is evident. During the LGM and the HS1, the contribution to evapotranspiration was maintained at a high level─ca. 80%, while only 20% of precipitation was converted to runoff. With climate warming, the contribution decreased to ca. 55% for evapotranspiration and increased to ca. 45% for runoff from the BA to the EH. Then, the contribution to runoff returned to ca.35% in the LH.

Fig. 4. Cycles in the different periods. The two detected cycles with most and secondary significant power spectrum density (PSD) are marked red and green, respectively.

Fig. 5. ISM variation for the YZRB over the past 22 ka. Indicators include simulated meridional wind speed anomaly, integrated ISM proxy based on four cave δ18O records Sahiya cave in Northern India(5.7-0 ka BP), Hoti cave in Northern Oman (9.6-6.0 ka BP), Mawmluh cave in Northeastern India(33.8-5.5 ka BP), Qunf cave in Southern Oman (10.6-2.7 and 1.3-0.3 ka BP) (Liu et al., 2020b), and the δ18O series for 5-22 ka BP from Indian Mawmluh Cave (Dutt et al., 2015).

3.4. Deglaciation related to abrupt runoff change

The glacier area has decreased from the LGM to the present (Fig. 8). The most significant decline occurred in the latter half of HS1. The glacier area ratio for the whole basin decreased from ca. 3 to ca. 1.5. Retreating glaciers during the late HS1, which converted to runoff, correspond to the fast increase of runoff (Fig. 2). However, subregions of the YZRB behaved differently. The glacier area in the middle reaches, which decreased by ca. 60%, shows the most evident change. Retreat rates of glacier area are ca.50% and ca. 30%, for the upper reaches and the lower reaches, respectively. The YZRB was subject to the YD cooling event; the response of glacier area is still smaller than that during the HS1. From the BA to the EH, there is a second significant decline in ice cover for all the subregions. For some extreme river runoff events, when surface temperature increased, warming could reduce the dam stability and glacier-meltwater changed the river runoff to flow into a dammed lake, suggesting that deglaciation, especially between the late HS1 and the BA, may have accelerated dam breaching to foster a megaflood.

Comparing ice cover with the simulated surface temperature, we found that the glacier area degrades generally as the temperature increases (Fig. 8). However, the direct correlation between these two parameters is not evident. When cooling occurs, there is no increase in the glacier area, especially during the rapid cooling periods, such as the YD and the early HS1. During the YD period, the fast drop in surface temperature can reflect climate change but it was not reflected in reducing the deglaciation.

Other than the surface temperature, the precipitation determining the water source of a glacier is also an important factor for glacier change. The standardized partial regression coefficients were used to show the relative effects of precipitation and temperature change on the glacier area(Fig. 9). For the whole basin, the effects of precipitation are very close to those of temperature during all the periods(Fig. 9a). After the YD cooling event, there are mainly positive precipitation anomalies, and the amounts of precipitation remain stable (Fig. 2). With certain continuous water recharge, the temperature gradually commands glacier change. Nonetheless, the situation in different regions of the YZRB are not the same (Figs. 9b-d). At the LGM, the glacier change in the upper reaches was dominated by temperature, while that in the lower reach was dominated by precipitation. Relative effects were stable in the HS1, indicating that the contribution of precipitation was slightly less. During the BA-YD period, only in the upstream region did the surface temperature affect the glaciers more than the precipitation. However, the result for the Holocene is contrary to that of the BA-YD period.

Fig. 6. Correlations between decadal-mean annual precipitation and runoff anomaly for the YZRB in different periods.

Fig. 7. Contribution of decadal-mean annual precipitation to evapotranspiration and runoff.

Fig. 8. Changes in glacier area and surface temperature. Glacier area ratios of transient glacier area to the current one are shown for glaciers in the upper, middle, lower reaches, and the whole river basin,respectively. Simulated surface temperature is used to reflect the thermal environment.

3.5. Discussion

The periods of unstable climates, were short but were capable of triggering major land-forming events such as megafloods and giant landslides (Clarke et al., 2003; Bookhagen et al., 2005; Kumar et al., 2014). Bookhagen et al. (2005)show that fivefold increases in sediment-transport rates recorded by sediments in landslide-dammed lakes characterized episodes of high climatic variability between 10-4 ka BP.Srivastava et al. (2017b) also found that most extreme floods induced by glacial lake outbursts occurred when the ambient climate was significantly wetter. The Hydro-climatic information of the current study may promote further understanding of high-impact water-related events in the YZRB.

Differences between the variations of runoff and the ISM proxy are significant during the YD (Figs. 2 and 5).The runoff of Tso Kar lake basin (Ladakh, India) decreased sharply from the BA to the YD as influenced by the monsoon (Wünnemann et al., 2010). The evolution of the ISM in the Bengal region also indicates a weak signal in the YD(Contreras-Rosales et al., 2014). Lake levels atop the Tibetan Plateau also were unusually low (Günther et al.,2015; Hou et al., 2017). However, the cooling effect on total runoff of the Yarlung Zangbo River was not detectable in our reconstructed data. Furthermore, the proxybased weaker ISM after the LGM caused by the Heinrich event has not been identified by the climate model (Fig. 5).This outcome may be due to the unmatched spatial scale between the global model and the local climate at the cave sites, or other incomplete process simulations because the trends of ISM reconstructed from different locations are not coincident (Hong et al., 2003; Contreras-Rosales et al.,2014; Dutt et al., 2015). The other explanation may be that the amplitude of the simulated global climate cooling for some period via the Community Climate System Model is less than that reconstructed from Greenland ice cores (Liu et al., 2009).

Fig. 9. Standardized partial regression coefficients of precipitation and surface temperature contributing to glacier area change for (a) the whole basin, (b) the upper reaches, (c) the middle reaches, and (d) the lower reaches. The larger the absolute value of coefficient, the greater the effect on the dependent variable.

In terms of glacier dynamics, glaciers change non-linearly with precipitation and temperature, which are not only related to climate but influenced by the topography and the physical properties of the glaciers. In most cases, data on paleotopography are absent when implementing fine-resolution glacier modeling. Furthermore, the physical properties of glaciers in ungauged regions are determined empirically.Precipitation-runoff relation may be altered significantly due to the regulating effects of glacier behaviors. At the basic level, the regulation effects depend upon the contribution ratio of the glacier melt to river runoff. However, finescale information such as daily temperature and precipitation to estimate glacier melt at the basin scale is difficult to obtain, due to limited model resolution at present and lack of direct measurements in the field.

On-going plateau climate change and fluvial erosion in the YZRB destroy many geological indications, inhibiting the paleo-hydrological record reconstruction. Reconstructing reliable climate proxies in the YZRB, harmonizing existing proxies from different sources in surroundings of the basin, and implementing high-resolution paleoclimate modeling in future work will allow the spatially-explicit analysis of hydro-climatic characteristics and their relationships with water-related extreme events.

4. Conclusions

In this study, we analyze multi-source climatic and hydrological data to elucidate the response of glaciers and runoff in the YZRB to climate change since the LGM. The mean river runoff stayed around a low level of 145 billion cubic meters per year until an abrupt increase in the BA. The annual runoff reached a maximum of 250 billion cubic meters in the early Holocene and then decreased to the current value of 180 billion cubic meters at a rate of 6.4 million cubic meters per year. The low runoff during the LGM and HS1 is likely attributed to such a small contribution of precipitation to runoff and large glacier cover. The percentage of precipitation to runoff was only 20%. The fastest deglaciation occurred during the late HS1, when nearly 60% of glacier area disappeared in the middle reach, 50% in the upper reach, and 30% in the lower reach. Deglaciation and increasing runoff between the late HS1 and BA may accelerate icedam breaches and lead to outburst flood events.

Acknowledgements.This work was supported by the National Natural Science Foundation of China (Grant No.91747207) and the project of CAS “Light of the West”. We appreciate Prof. Qing YAN for sharing the simulated glacier areas, Prof.Jianbao LIU for sharing the integrated ISM proxy, and Prof. Anil K. GUPTA for sharing δ18O records. TraCE-21ka datasets can be downloaded from the Climate Data Gateway at https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm3.trace.html. All other data drawn in this paper can be requested from the first or corresponding author.