# **Application of Climatic Data in Hydrologic Models**

Edited by Mohammad Valipour and Sayed M. Bateni Printed Edition of the Special Issue Published in *Climate*

www.mdpi.com/journal/climate

# **Application of Climatic Data in Hydrologic Models**

# **Application of Climatic Data in Hydrologic Models**

Editors

**Mohammad Valipour Sayed M. Bateni**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Mohammad Valipour Metropolitan State University of Denver USA

Sayed M. Bateni University of Hawaii at Manoa USA

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Climate* (ISSN 2225-1154) (available at: https://www.mdpi.com/journal/climate/special issues/ climatic data hydrologic models).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5065-7 (Hbk) ISBN 978-3-0365-5066-4 (PDF)**

Cover image courtesy of Dr. Mohammad Valipour

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

# **Contents**


### *Editorial* **Global Surface Temperature: A New Insight**

**Mohammad Valipour 1, Sayed M. Bateni <sup>1</sup> and Changhyun Jun 2,\***


This paper belongs to our Special Issue "Application of Climate Data in Hydrologic Models". Here, we represent an example of the importance of climate data in terms of global surface temperature. This paper investigates the changes of 20-year (2000–2019) mean surface temperature (ST), wind speed (WS), and albedo (AL) data from the Global Land Data Assimilation System (GLDAS) over the globe with respect to those in 1961–1990. Moreover, we assess if the alterations in ST are affected by the changes in WS and AL. We also discuss the main reasons for the variations observed in ST, WS, and AL on global and hemispheric scales.

Our planet is within ~1 ◦C of its highest surface temperature in the past million years [1]. The most important indicator of global warming and climate change is the global mean surface temperature (GMST) [2]. Figure 1a shows the time series of monthly GMST anomalies from the Hadly Center/Climate Research Unit (HadCRUT) [2,3]. The mean of monthly GMST anomalies in 2000–2019 is 0.54 ◦C higher than that in 1961–1990 (Figure 1a) (see also [1,4]). The standard deviation (SD) of the decadal variation of GMST anomalies is reduced by 0.05 ◦C in 2000–2019 compared to 1850–1899 (Figure 1a). This is due to the climate sensitivity to oscillation indices, particularly El Niño/Southern Oscillation (ENSO) [5]. Indeed, there were more ENSO events in 1850–1900 than in 2000–2020, which affected the GMST anomalies. For example, the high positive GMST anomaly in 1878 (shown by the black circle) is caused by the late 1870s El Niño (Figure 1a) [1].

Figure 1b shows the variation of Global Land Data Assimilation System (GLDAS) ST data (ΔST) from 1961–1990 to 2000–2019 over the globe. The results show the GMST increase of 0.66 ◦C from 1961–1990 to 2000–2019. Analogously, Hansen et al. [1] reported the GMST growth of ~0.2 ◦C per decade during the last three decades. As indicated in Figure 1b, in general, the Northern Hemisphere (NH) has a higher variation of ST (ΔST) than the Southern Hemisphere (SH) [4]. The averages of ΔST in the NH and SH are 1.27 ◦C and −0.18 ◦C, respectively. Increasing greenhouse gas (GHG) emissions and variations of the North Atlantic Oscillation (NAO) are the main causes of increasing ST across the globe and particularly in the NH [1,6]. Effective policies must focus on using clean energies to reduce the GHG emissions. Such policies can facilitate achieving the 2015 Paris COP21 Agreement to maintain the GMST below 1.5–2 ◦C above pre-industrial levels. The decreasing ST in the SH is associated with the variations of the SH Annular Mode (SAM) and intensification of the South Pacific Anticyclone [7,8]. Fogt and Marshall [7] showed that adiabatic cooling over the Antarctica due to the variations of SAM aids in decreasing ST in the SH.

Figure 1c indicates the boxplot of GMST anomalies in each season. The horizontal line within the box indicates the median. The upper and lower edges of the box represent the 75% and 25% percentiles, respectively. The upper and lower ends of the whiskers show the maximum and minimum values, respectively. Outliers are observations beyond the end of the whiskers. As shown, the variation of the GMST anomaly during December–February (DJF) is more than that of other seasons. Vose et al. [9] also observed the highest upward

**Citation:** Valipour, M.; Bateni, S.M.; Jun, C. Global Surface Temperature: A New Insight. *Climate* **2021**, *9*, 81. https://doi.org/10.3390/cli9050081

Received: 27 April 2021 Accepted: 7 May 2021 Published: 12 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

decadal trend of ST in the NH (particularly high latitudes) in DJF when the minimum ST often occurs. In addition, the global minimum ST has increased approximately twice as fast as the global maximum ST since 1950 [9]. Hence, a larger variation in the GMST anomaly is seen in DJF (Figure 1c).

Comparison of Figure 1b,d illustrates a larger variation of ST during DJF in both hemispheres. During DJF, the mean of ΔST in the NH is 1.8 ◦C (Figure 1d). However, in some parts of Eastern Russia, it can reach up to 20 ◦C (Figure 1d). Hansen et al. [1] concluded that the global warming of more than ~1 ◦C could lead to a dangerous climate change and impacts the sea water level.

From 1961–1990 to 2000–2019, the GLDAS ST (WS) increased (decreased) in the northwest of North America (Alaska) and northeast of Asia (Russia) (Figure 1b,e). Both ST and WS increased in Canada (except in the western part), most regions of Australia and North Africa. In most regions of Africa and South America, WS enhanced, but ST reduced. Figure 1f shows the changes in the ST-WS correlation (Δ*ρST*−*WS*) from 1961–1990 to 2000–2019. High positive or negative values of Δ*ρST*−*WS* (>1 or <−1) in Southern Africa, North Greenland, Western America, Eastern Russia, and Central and Eastern Asia denote the changes of vegetation cover [10].

As can be seen in Figure 1e, winds slow down in most of the coastal regions as well as in Central and Eastern of Asia due to an increase in the surface roughness. These findings are consistent with those of Vautard et al. [10] who reported the highest downward trend of WS in Central and Eastern Asia from 1979 to 2008 due to the land-use change and/or biomass increase. They demonstrated that increasing the Normalized Difference Vegetation Index (NDVI) in Eurasia could explain 25–60% of the NH atmospheric stilling.

The increase of AL plays an important role in future climate change [11]. Ghimire et al. [12] reported an increase of about 1% in the global albedo due to the land cover change during 1700–2005. That caused the global atmospheric radiation to vary by almost 0.15 W/m2. Figure 1g represents changes in the GLDAS AL from 1961–1999 to 2000–2019 over the globe. On average, AL increases by 3.19% in the NH. AL varies with changes in cloud fractional coverage, aerosol amount, and land cover type [13]. If the variation of AL is caused by changes in aerosols and land cover type, ST reduces by rising AL [13]. In the NH, AL increased from 1961–1999 to 2000–2019 (Figure 1g), while ST did not decrease (Figure 1b). Therefore, the increase of AL is indicative of climatologically significant cloud-induced variations in the Earth's radiation budget, which is in line with the GMST growth in recent decades [11]. Another reason of AL growth may be the intensified wildfires in the NH boreal forests (due to droughts), especially in Canada and Alaska (see also [14,15]). In this region, the ST increase (Figure 1b,d) along with the precipitation reduction led to drought, which enhanced the wildfire severity [14].

Through changing ST, NDVI, AL and evapotranspiration, wildfire is the dominant driver of water, energy, and carbon cycles as well as vegetation dynamics in the North American boreal region [15,16]. The aforementioned region has the highest increase of AL that can reach up to 36% (Figure 1g). Likewise, Jin et al. [16] indicated that the increased severity of wildfires in the North American boreal forests during 2000–2009 enhanced AL by ~60%. Potter et al. [15] showed that climate change counteracts the cooling impact of postfire AL growth in the NH boreal forests. This justifies the increase of ST despite rising AL in the NH (Figure 1b,g).

Figure 1h compares the changes of ST, WS, and AL from 1961–1990 to 2000–2009 in the NH and SH. Both ST and AL increased in 86% of the NH and 9% of SH regions. Increase of ST in most regions on the NH leads to drought events, which intensify wildfires, ultimately rising AL. In addition, variations of ST are consistent with those of WS (i.e., both increase or decrease) in 53% and 51% of the areas of NH and SH, respectively. This happens because WS can affect ST through modifying evapotranspiration [17,18].

Regional studies are necessary to improve our understanding of the variations in ST, WS, and AL [4]. Some important regions for further investigations are the High Mountains

of Asia, which demonstrate extreme values of ΔST (Figure 1b), ΔWS (Figure 1e), and ΔAL (Figure 1g) (see also [10,19]).

**Figure 1.** (**a**); Time series of monthly and decadal GMST anomalies. The 1850–1899 and 2000–2019 periods are shown by pink colors to compare the variations before the 20th century and after the 21st century, (**b**); variation of the GMST, (**c)**; boxplots of the GMST anomalies. The horizontal line within the box indicates the median (50% percentile). The upper and lower edges of the box represent the 75% and 25% percentile, respectively. The upper and lower ends of the whiskers represent the maximum and minimum values. Outliers are observations beyond the end of the whiskers, (**d**); variation of the GMST in DJF, (**e**); variation of WS, (**f**); variation of the correlation between the GMST and WS, (**g**); variation of AL, (**h**); Hemispheric variation of the GMST, AL and WS. (**a**,**c**) are plotted based on the HadCRUT-version 4.6 product [20]. (**b**,**d**–**h**) are drawn based on the GLDAS (1◦ × 1◦ spatial resolution and monthly temporal resolution) dataset [21].

At the end, we would like to thank all the authors for their contributions to this Special Issue. We hope that the papers published in this Special Issue improve our knowledge in terms of application of climate data in hydrologic models.

**Author Contributions:** Conceptualization, M.V.; investigation, M.V., S.M.B., C.J.; writing—original draft preparation, M.V.; writing—review and editing, S.M.B., C.J., supervision, S.M.B., C.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding authors upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Traditional Water Governance Practices for Flood Mitigation in Ancient Sri Lanka**

**Vindya Hewawasam 1,\* and Kenichi Matsui <sup>2</sup>**


**Abstract:** The tank cascade system, which emerged as early as the fifth century BC in Sri Lanka's dry zone, has been portrayed as one of the oldest water management practices in the world. However, its important function as flood management has not yet been thoroughly examined. In this paper, we argue that the main principle behind the tank cascade system is not only to recycle and reuse water resources by taking advantage of natural landscapes but also to control floods. This paper examines the evolution of traditional water management and flood mitigation techniques that flourished in pre-colonial Sri Lanka. This historical examination also sheds light on recent policies that exhibited renewed interests in revitalizing some aspects of the tank cascade system in Sri Lanka's dry zone. This paper shows how ancient Sinhalese engineers and leaders incorporated traditional scientific and engineering knowledge into flood mitigation by engendering a series of innovations for land use planning, embankment designs, and water storage technologies. It also discusses how this system was governed by both kingdoms and local communities. Water management and flood control were among the highest priorities in urban planning and management. The paper thus discusses how, for centuries, local communities successfully sustained the tank cascade system through localized governance, which recent revitalized traditional water management projects often lack.

**Keywords:** tank cascade system; dry zone; water governance; flood control; traditional knowledge; community participation; Sri Lanka

#### **1. Introduction**

Water management is one of the fundamental requirements for the survival and prosperity of civilizations. Historically each civilization developed its unique water management practices that reflected the surrounding topography, climate, soil and utility purposes [1,2]. The traditional water management system that developed in the dry zone of Sri Lanka more than 2500 years ago is one of the oldest known water management systems in the world [3]. This ancient hydraulic civilization uniquely engineered storage dams [4] and water distribution systems [5,6].

According to the *Mahavamsa*, an earliest known chronicle that depicted ancient Sri Lanka, when Indian Prince Wijaya landed Sri Lanka from India in 543 BC, he observed irrigation practices among Indigenous Sinhalese people. At the port of Mannar, where he landed, the prince found many water tanks (or reservoirs) with cool water that replenished a great garden [5].

In the twenty-first century, this traditional water management is still in practice to some extent although much of it has been disrespected due to the introduction of Western water management systems under the colonial regime as well as in the age of the more contemporary international cooperation regime. However, somewhat reversing this trend, the Sri Lankan government has recently acknowledged the importance of historically

**Citation:** Hewawasam, V.; Matsui, K. Traditional Water Governance Practices for Flood Mitigation in Ancient Sri Lanka. *Climate* **2022**, *10*, 69. https://doi.org/10.3390/ cli10050069

Academic Editors: Mohammad Valipour and Sayed M. Bateni

Received: 8 April 2022 Accepted: 11 May 2022 Published: 13 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

practiced water management and reinstalled some in rural regions for flood protection and irrigation purposes.

Some scholars have reexamined Sri Lanka's traditional water management system partly for the purpose of enhancing modern day climate resilience actions. Some emphasized the implications for drought mitigation and rain water harvesting [7–11], agricultural developments [5,12,13], ecosystem management [2,14–17], and socio-economic development [4,18–20]. Literature shows that a similar tank system existed in semiarid southern India, but its main purpose was to provide water for paddy cultivation [21,22]. Qanat is another traditional irrigation water management practices that existed in semi-arid regions of Morocco, Spain, Syria, Iran, and Central and Eastern Asia. Ancient societies developed underground networks for the transportation of water [1].

Looking at this growing trend of studies, what is missing is a linkage between traditional water management and flood mitigation practices in Sri Lanka and elsewhere. Some researchers did mention about traditional flood mitigation functions in Sri Lanka [2,13,22–25], but the question remains as to the extent to which traditional water management practices were systemically arranged for improving or supplementing flood protection. In other words, we argue, flood mitigation has long been integral part of Sri Lanka's water management system.

This paper, therefore, seeks to understand Sri Lanka's traditional flood mitigation functions and technologies that evolved through time in its dry zone. This said, some may argue that this type of examination requires hydrological modeling or engineering investigation to truly understand the effectiveness of ancient flood mitigation infrastructure [26]. However, our main focus in this paper is rather to trace how past practices took shape in time, given urgent needs of local Sri Lankan farmers to mitigate flood risks. The IPCC report and other recent studies on climate change adaptation and disaster mitigation emphasized the importance of better understanding locally developed adaptation and mitigation practices as a way to enhance local disaster response capacity and participation [6,13,22,27–30].

In the following discussion, we first look at the development process and functions of the traditional system. Then, we examine water and flood management practices in ancient cities. Finally, we discuss the water governance system and its sustainability. For our examination, we used historical records, secondary sources, institutional reports, and audiovisual sources. In March 2019, the authors visited several ancient water management sites, including Sigiriya and Polonnaruwa to collect documentary and visual information. We also collected information in Colombo in the same year to find out what has already been known in the country about its ancient water management system.

#### **2. Development of the Tank Cascade System in the Dry Zone of Sri Lanka**

In Sri Lanka's dry zone annual mean precipitation is about 1750 mm whereas annual mean evaporation ranges from 1700 mm to 1900 mm [31]. About 80% of the annual rainfall occurs during the northeast monsoon season from November to February when flash floods often occur. Seasonal rivers and so-called *Villu* (wetland ecosystem in floodplains) are natural water bodies that emerge during these months. The earliest inhabitants were recorded in the lowland areas such as Anuradhapura and north central parts of Sri Lanka in the ninth century BC [32,33]. They lived along rivers and water bodies, collecting and storing water partly for drinking and irrigation purposes [2,6].

The earliest available information on hydraulic civilization of Sri Lanka dated back to the sixth to fifth century BC [4,5]. Archaeological studies show a network of tanks that were interconnected by streams and waterways. Today, this water network is commonly known as the tank cascade system [5,8,9,34–36]. The main hydrological principle behind the tank cascade system is to recycle and reuse water through a network of small to large scale tanks within a catchment. It also considers storing, transporting and distributing water for mitigating floods and droughts [2]. The International Union for Conservation of Nature in Sri Lanka identified four main functions of these ancient structures: (1) capturing rainwater to minimize floods; (2) storing rainwater; (3) recycling used water; and (4) mitigating drought impact [25,27]. Other than these, the cascade system sustained the local ecosystem as ancient engineers carefully used natural landscapes to enhance water storage [37].

The flood mitigation of the tank cascade system entailed engineering techniques of water and sedimentation flow control along with the protection of banks from erosion. Ancient Sinhalese people constructed granite structures and pillars. In order to protect the embankment from breaching and flooding, with improved technologies for metallurgy, iron was used to strengthen the structure [5]. The knowledge of iron metallurgy was introduced to Sri Lanka as early as the tenth century BC. The archaeological sites of Aligala, Sigiriya and Anuradhapura show evidence of iron smelting in the ninth century BC [38]. The construction of large tanks emerged in the fourth century BC. At the time of increasing water levels, tanks had spillways to safely release excess water from one tank to another. Small tanks were built in low plains between hills by connecting them with embankments [36]. For example, *Basawakkulama* tank, the earliest recorded large-scale tank was built in the Malwathu Oya Basin in about 430 BC. With over 3000 feet in length and 21 feet in height, the dam had storage capacity for cultivating 350 acres. A large number of small tanks were built in the same basin to avoid possible disasters from flooding [4].

Mahatantila et al. [39] identified three main components of tanks: (1) upper periphery, (2) bund/embankment and (3) tank body. However, in the following discussion, we add one more component, which is especially important for flood management; that is, the lower periphery of the tank where human settlements with paddies were located. Figure 1 shows how these components were typically laid out. Paddy cultivation was the main livelihood practiced by early inhabitants. The paddy cultivation of Sri Lanka dates back to the ninth to sixth centuries BC. During this period, ancient farmers domesticated cattle. Cattle were used for harrowing paddy fields [32,33,40].

**Figure 1.** Main components of the tank.

In the upper periphery of the tank running water in streams was filtered through forests with patchy water holes or bogs. Rain-fed farms were located here [2]. The ancient law prohibited felling trees as forests were important to manage water quantity and quality. Ancient people developed water holes (*godawala*) partly to prevent sediments from entering the tank [19,23]. Below these water holes, a water filtering area called *perahana* was created with water grasses like reeds [23,27].

Ancient Sinhalese protected the embankment from wind, heavy rain and waves by building stone liners on embankment walls [2,20]. The tank embankment was basically made of earth and granite rocks. Large-scale embankments with 30–40 feet deep reservoirs consisted of unique and intricate engineering innovations. The height of the embankment was carefully designed not to flood the upper stream area [14]. The embankment was installed with a sluice gate (*sorowwa*), valve pit (*bisokotuwa*), water level indicator (*diyakata pahana*), spillway (*pitawana*) and embankment protector (*ralapanawa*). The main purpose of sluice gate (*sorowwa*) was to regulate water release without flooding lower stream areas. The water level indicator helped decide when to release water. The valve pit or *bisokotuwa*, which was attached to the sluice gate on the bottom of the embankment, was basically a rectangular buffer room that was created to temporarily gather water from the lake through the sluice. When the level of gathered water in the room was raised above the sluice gate level but below the reservoir water level, water was released toward the lower periphery through the other gate(s). The location and arrangement of these water gates differed by region and embankment, showing engineering diversity [16,41]. The *bisokotuwa* is still in operation at Kalawewa, one of the largest tanks built in 477AD during the period of Anuradhapura kingdom [2].

In the upper edge of tanks, a so-called tree belt (*gasgommana*) had a number of planted trees partly to protect the embankment. It also provided the habitat for fish and other aquatic species [2,23]. The trees became partly submerged in water during the heavy rain period. This tree belt also acted as a wind barrier and reduced the waves in the tank.

In the lower periphery of the tank, when water was released from a sluice gate through a valve pit it ran through an interceptor (*kattakaduwa*). The interceptor is a reserved land for the purpose of controlling soil erosion and water contamination. Villagers took drinking water below the interceptor. It also provided water to farms. The surrounding village was protected from water inundation with a hamlet buffer area (*thisbambe*), shrub land (*landa*) and drainage (*kivul ela*). The trees that have high heavy metal and salt absorption capacity with a strong root system were planted [2,23]. Being in the high elevated areas near the interceptor, villagers could observe flooding or damages to the embankment. The hamlet was surrounded by the hamlet buffer area that was used for common perennial cultivation (e.g., mango, coconut) and resting places for buffaloes [20,23]. Paddy fields below the interceptor functioned as wetlands during heavy rains to keep temporary flood water. When water is not enough for the whole paddy lands, all farmers cultivated equally (*Bethma cultivation),* limiting the paddy area to be irrigated [13,15]. The villagers used the shrub land for home gardening, such as chena cultivation. The excess water of the paddy fields flowed to the drainage area that was used for common village purposes to absorb salt and other contaminations [23,27]. Through the drainage and other natural streams, water reached the next tank.

#### **3. Flood and Water Management Techniques in Ancient Kingdoms in Sri Lanka**

In planning cities, villages and monastery complexes, ancient Sinhalese engineers carefully considered water sources, water uses and landscapes [42]. Flood control is one of the main requirements of city planning. Figure 2 shows the locations of the ancient kingdoms, main tanks and rivers. The historical records show that water was used not only for drinking and irrigation purposes but also for public bath and recreational activities. Traditional knowledge on rainfall patterns, land use planning and landscape helped ancient people maximize the use of water resources [17].

**Figure 2.** The locations of the ancient kingdoms, main tanks and rivers.

For example, in the early fifth century BC, the city of Anuradhapura adopted the ancient Indian "Nandyavarta plan" for water resource management [42]. It laid out the city in a circular shape. The central area was surrounded by walls with four gates that faced each cardinal direction. Anuradhapura was established along the Malwathu Oya (or river), the main river of the area. Anuradhapura's engineers constructed five large tanks around the city in about 437 BC. In order to protect the city from floods and droughts, they also constructed many small tanks in the same valley. Between the river and the tanks, three green parks were established mainly for recreational purposes [42]. These parks acted as water retention facilities during floods. Villages and king's palace were located below *Basawakkulama* tank. Its L-shaped embarkment was designed to take advantage of the surrounding landscape. It supplied both drinking and irrigation water. Bathing also became important part of Anuradhapura residents' lifestyle [42].

Sigiriya fortress, the present-day UNESCO World Heritage site, is another example to show a complex outlook of water management in an ancient city. The annual rainfall is the only possible water source here [8,17]. The fortress, which was built on the top of a gigantic rock as well as its surrounding areas, was built by King Kassapa in the fifth century AD (477–495). Human settlements began in this area as early as the third to the first century BC [43–45]. These early settlements were basically for monks who lived in caves of Sigiriya Rock. The rock walls just above caves were carved out like a gutter to keep out rainwater from flooding dwelling areas [43,45]. The cave entrances were then plastered for further protection. Archeologists identified about thirty such locations [43].

Later, King Kassapa developed an urban complex here [43,45]. Residents took water from the Sigiriya Oya and stored in a tank near Sigiriya Rock. Engineers at the time built storage tanks, cisterns, water-courses, underground and surface drainage to managed water in the city. All storage ponds and bathing pools were paved with marbles and pebbles to enhance water retention. In addition, natural depression areas were used to collect rainwater. The city was designed to control flow velocity, runoff discharge, and flow distance. For example, non-structural depression areas and drainage patterns were used to direct rainwater to the structural ponds located in lower elevations. During the process, water was filtered and velocity was reduced to control soil erosion [17]. Sigiriya also had water gardens and fountains [17,45]. The pools in the water gardens and Sigiriya tank were interconnected through underground drainage. This helped fill the pools automatically [8,42]. During the rain, the water garden can function as a water storage facility. In maintaining pools for bathing purposes, water was supplied from storage tanks, and the used water was released to moats through a separate drainage. The fountains were connected to special underground channels [42]. The moats were located in the lowest area of the land and excess water flows into the moats by reducing floods (Figure 3).

**Figure 3.** A water storage tank in the Sigiriya complex, Photo courtesy: Kenichi Matsui.

King's palace and the rock garden were located on the three-acre rock summit area about 360 m above the sea level [43]. Roofs were designed to collect and transport rainwater to the main water storage area. A drainage outlet was constructed to dispose excess water and prevent flooding. The surface area was terraced with the western side as the highest point [42].

After the demise of Sigiriya city, the Anuradhapura kingdom was reestablished in the 5th century AD about 80 km northwest from the rock. In its monastery site called Abayagiri, twin ponds were created in a low-lying area partly for monks to bathe. It is considered one of the best hydrologic engineering marvels of ancient Sri Lanka [46]. Underground pipelines were established to connect the ponds to Tissawewa, Basawakkulama, and Nuwarawewa tanks around the city by drawing water from the Malwathu Oya [42]. These pipelines ran through a number of small sediment/debris control tanks [47]. An enclosing wall was built around the ponds to control the possible spillage [46]. Also, wastewater outflows ran through wetlands for purification. Then the water was released back to the same river [8]. Each component of the ponds was carefully designed to protect the monastery complex from flooding.

Wastewater is a significant threat to health, particularly during flood events [28,48,49]. Ancient Sinhalese developed and practiced wastewater management. In Anuradhapura and Polonnaruwa different types of lavatories were developed [50]. Here urinals were collected in pits through terracotta pipes. Sands, lime powder and charcoals were used to purify wastewater [8,50]. In some places, separate septic tanks were used to store

wastewater. These places were kept some distance from residential areas in order to avoid wastewater spillage during floods.

In the reign of King Parakramabahu (1153–1183 AD), water management technologies were further refined [5]. His engineers built several large-scale tanks and irrigation systems [23]. They constructed more than 163 major tanks, 2376 minor tanks, 165 anicuts and 3910 diversion channels [51]. The capital was located in present-day Polonnaruwa city, about 50 km east from Sigiriya. King's palace was located very close to the tank called *Parakyama Samudraya* (Figure 4) that had nine-mile-long embankment. In order to protect the palace from floods, huge brick walls were constructed along the tank side of the palace. There were several non-structural natural drainage facilities inside the walls to temporarily store water. The ground was also covered with grass to reduce soil erosion and trap debris. A few sluice gates were installed along the brick walls. A drainage canal was installed in the other side of the palace. *Kumara Pokuna* near King Parakramabahu's palace was one of bathing ponds that might have functioned as flood control structure (Figure 5). Similar to the Anuradhapura pipeline system it was connected with several drainages to purify water [17]. Even today this system is functioning well.

**Figure 4.** View of *Parakrama Samudraya* Photo courtesy: Kenichi Matsui.

**Figure 5.** *Kumara Pokuna* in Polonnaruwa Photo courtesy: Kenichi Matsui.

#### **4. Water Governance in Ancient Sri Lanka**

The ancient chronicles of *Mahavamsa*, *Dipavamsa* and *Culavamsa* as well as remaining cave, rock and slab inscriptions show evidence of Sri Lanka's water governance from the fourth century BC to the thirteenth century AD [3]. Some of these documents tell us how early Anuradhapura kingdom governance practices emerged with professionals and water ownership. The government imposed an income tax and other rules on using water [52]. These rules basically relied on community participation and involvement [10,13,17]. Later, the governance system gradually changed from a community-based system to a centralized one although the small-scale tank cascade system remained under community management [3].

From the fourth century BC to third century AD, water rights in general were held by individuals, kings, elites, local chiefs and families [3]. Kings and elites could grant water rights to monks mainly in the period between the second century BC and eighth century AD. Buddhist temples then administered water allocations. Until the second century BC, the management of tanks, including flood control and maintenance, were mainly undertaken at the village level, which mainly consisted of farmers [3,13]. Works for repair, desiltation and cleaning of the tank during the dry season were shared among farmers proportionately to the land ownership. Each farmer provided his or her service free on certain days [20]. The community also planted trees to strengthen the stability of the tank embankment and the interceptor (*Kattakaduwa*). In times of water shortage, the *Bethma* rules made farmers share water for paddy cultivation. This set of rules are still practiced today among some farmers during water shortage [13]. Here the village head or prominent leader decide the area for cultivation each year based on water availability in the tank. Farmers then received equity-based water allocation based on land ownership [13,20,40,53].

In the second century BC, localized water governance was gradually replaced with a centralized system. Different professions emerged as a result, such as flow operators (*Naguli*), canal officials, and proprietors of ferries (*Parumaka Thota Bojhaka*), and proprietor of tanks (*Parumaka vapihamika*). This institutionalized governance made it possible to sustain the food supply of a large population [3]. In the ninth century AD, the last phase of kingdom of Anuradhapura further institutionalized water governance by establishing specialized committees to maintain large-scale tanks [3].

Water rites, rituals, and customs played an important role in ancient water governance. For example, the king granted water rights through the "ceremony of golden vase", in which water was poured from a golden vase into farmers' hands [3]. The king also participated in festivals that sent a signal to commence ploughing, sowing, and harvesting. *Pen Pidima* ceremony offered fresh water of a tank to Buddha statues to pray for fertility [18].

#### **5. Abandonment of Dry Zone Water Governance to Contemporary Water Governance**

By the mid-13th century AD with the collapse of the Polonnaruwa kingdom, the centralized large-scale tank cascade governance was largely abandoned in many parts of the dry zone [3,4,31,54]. Although community-driven small-scale cascade systems remained in practice with varying degrees until the end of the 15th century [31]. European colonization under the Portuguese, Dutch, and British from the 16th through 18th centuries systemically and gradually disempowered traditional local authorities for water governance [4].

In 1832, about 17 years after the British established its colonial government in Kandy, it abolished the Kandyan *rajakariya* system, which imposed compulsory labor for public works, claiming that it resembled a form of slavery [3,31]. At the time many villages governed local affairs through *Gansabhawa*, a council composed of representatives of villagers. This council depended largely on the availability of village labor under the *rajakariya* system. As the British colonial regime further tightened restrictions on it, canals and other traditional water management works gradually fell into serious decay [54]. This change led to the deterioration of the community tank cascade system in many parts of the dry zone [3,31].

The British then empowered the Temple Land Commission (1856) and the Service Tenure Commission (1870) to govern tanks and surrounding villages where village heads used to control local affairs [3]. From 1870 to 1897, the colonial government repaired and restored a large number of village tanks, including *Kalawewa*, one of the largest tanks built in 477AD and Giant's Canal (*Yoda ela*), which was built in 459AD. Village communities provided their labor in voluntary basis for these tasks [54]. In 1900, as more educated Sri Lankan elites had been brought into civil services by then, the colonial government hired many of them and established the Department of Irrigation, which took control over public works. This department remains as one of the oldest departments in Sri Lanka. One of its mandates was flood protection [4,55]. During its 50-year operation, it restored almost all large-scale tanks and anicuts [4].

Soon after independence in 1948, the government of Sri Lanka attempted to improve water management. For example, in 1952, it constructed the nation's largest reservoir called *Senanayaka Samudraya* [29]. As the Gal Oya basin often experienced floods and droughts [55], disaster mitigation was one of the main components of the project. However, the project ultimately failed to achieve its initial goals due largely to top-down decision making where community participation was not incorporated [29].

In the 1950s, the Sri Lankan government began to place more emphasis on regional tank cascade governance. The Paddy Lands Act of 1958, for example, authorized the Department of Agrarian Services to maintain all small-scale water works. This act placed village tank cascade systems directly under the Department of Agrarian Services, a central government authority. It led the restoration and rehabilitation of small-scale village tanks including canals and flood mitigation structures [31].

In 1977, the Sri Lankan government undertook a comprehensive basin development called the Mahaweli Accelerated Development Project in the Mahaweli River Basin and created the Mahaweli Authority [4]. Flood control in the lower Mahaweli River was one of the main objectives [56,57]. It constructed new Western-style water reservoirs such as Kotmale, Victoria, and Maduru Oya. Kotmale reservoir was specifically designed for flood control by making it possible to transfer flood water from Polgolla to the dry zone where restored ancient tanks are located (e.g., *Parakrama Samudraya*, *Minneriya*, *Kantale*, *Kaudulla*, *Kalawewa and Giritale*) [58]. After the construction of this reservoir was completed in the mid-1980s, the flood inundation risk in the downstream of the Mahaweli River was significantly reduced [58]. The Mahaweli Authority also reestablished ancient river connections among Sudu Ganga, Amban Ganga, and Dambulla Oya in central province [59].

In the 1980s, the government undertook projects to restore ancient water management by mobilizing local people. In 1981, the Gal Oya Left Bank Rehabilitation Project hired community labor and collected local knowledge about water management [29,60]. It created farming organizations to control local water use for domestic and agricultural purposes. It also funded channel maintenance by these organizations. As farmers in the dry zone are the ones to experience flood damage to their crops, these farmer organizations were expected to play active roles in flood management [29]. In 1982, for example, the Village Irrigation Rehabilitation Project and National Irrigation Rehabilitation Project aimed to repair and maintain minor tanks with community participation [37].

In 1991, the Agrarian Services Act induced the concept of joint water management between farmers' organizations and a government agency [31]. In order to promote participatory planning and management, stakeholders were engaged in *kanna* meetings (pre-seasonal meetings of farmers). Even today, these meetings are the most important decision-making bodies in operating local tanks. Their tasks include the joint maintenance of flood control bunds, sluices, and channels [31].

Post-independent tank rehabilitation programs were mainly for repair, maintenance and physical development of individual tanks rather than the whole network of the tank cascade system. This shortcoming resulted in tank sedimentation, water leakage, land degradation, biodiversity loss, and floods and soil erosion [37]. In the 1990s, the Shared Control of Natural Resources in Watersheds Project adopted an ecosystem approach through

community-based participatory watershed management. It expressed its strong interests in environmental sustainability and productivity improvement [60]. The project was implemented in northcentral province and southern province. It introduced farming companies to work independently in watershed management [61]. The project took the novel approach to address watershed issues by creating mini-projects among local communities and NGOs. This provided a sense of ownership for the parties involved in the project [62]. For example, in the Nilwala watershed, deforestation in the upper stream often flooded the lower basin due to sedimentation. To prevent soil erosion and sedimentation, the project attempted to protect existing forests and rehabilitate degraded forests. At villages, agroforestry practices were encouraged [63].

#### **6. Conclusions**

This paper has examined how Sri Lanka's traditional flood and water management evolved in its dry zone. Our goal was not to suggest that all traditional forms of flood control and community participation were effective. However, to better understand how Sri Lanka's flood management practices exist today in a hybrid form, which has incorporated traditional, colonial and modern technologies and practices, we found it imperative to clarify historical changes in flood management practices. Without knowing this complexity, locally viable flood mitigation measures cannot be sustainable. Therefore, it is important to understand how, in Sri Lanka, societies developed unique flood control practices, including governance frameworks and laws to ensure safety from disasters and equitable access to water. The tank cascade system is the result of their traditional knowledge on watershed and disaster management.

To demonstrate some recorded practices, this paper looked at how the tank cascade system captured rainwater, stored it, minimized flood impacts, maintained public health and conserved/nourished biodiversity. It somewhat resembled modern-day integrated water resources management. The traditional knowledge regarding engineering and metallurgy evolved and resulted in large-scale embankments for flood control along with sluice gate (*Sorowwa*), valve pit (*Bisokotuwa*), water level indicator (*Diyakata pahana*), spillway (*Pitawana*), and embankment protector (*Ralapanawa*).

Anuradhapura, Sigiriya, and Polonnaruwa cities developed both structural and nonstructural water infrastructures. Multiple-purpose structures were designed for flood mitigation, irrigation, purification, drinking, and recreation. For centuries, the tank cascade system was largely governed by the community. Experienced community leaders played a vital role in decision making. Community voluntary support for managing the commons is an important feature in water governance. Rights to water resources were shared among elite groups, monks, community and individuals. The development of water professionals, taxes and rules in managing water systems emerged when centralized water governance under kingdoms began to exert more control over water resources. Along with these institutional development, monks and people developed water rituals and customs that considerably influenced traditional water governance.

European colonial regimes, however, gradually eroded this intricate water governance practices. The abolition of the *rajakariya system* under the British rule led to the disuse of communal tanks as it became difficult to obtain local labor. After independence, the government of Sri Lanka showed its renewed interests in traditional water governance and undertook several large-scale water management projects, such as the Gal Oya Irrigation Project and Mahaweli Accelerated Development Project with renewed interests in locally viable traditional water management. In the 1990s, Sri Lanka's watershed restoration policy began to emphasize community participation and led to some positive results. More water governance projects are planned to take advantage of traditional systems and mobilize local participation although the overall impact of their effectiveness for flood control under escalating climate change conditions remain to be determined in the future.

**Author Contributions:** Conceptualization, V.H. and K.M.; methodology, V.H.; software, V.H.; validation, V.H. and K.M.; formal analysis, K.M.; investigation, V.H. and K.M.; resources, V.H. and K.M.; data curation, V.H. and K.M.; writing—original draft preparation, V.H.; writing—review and editing, K.M.; visualization, V.H.; supervision, K.M.; project administration, K.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Climate Change Impacts on Groundwater Recharge in Cold and Humid Climates: Controlling Processes and Thresholds**

**Emmanuel Dubois 1,2,3,\*, Marie Larocque 1,2,3, Sylvain Gagné 1,2 and Marco Braun <sup>4</sup>**


**Abstract:** Long-term changes in precipitation and temperature indirectly impact aquifers through groundwater recharge (GWR). Although estimates of future GWR are needed for water resource management, they are uncertain in cold and humid climates due to the wide range in possible future climatic conditions. This work aims to (1) simulate the impacts of climate change on regional GWR for a cold and humid climate and (2) identify precipitation and temperature changes leading to significant long-term changes in GWR. Spatially distributed GWR is simulated in a case study for the southern Province of Quebec (Canada, 36,000 km2) using a water budget model. Climate scenarios from global climate models indicate warming temperatures and wetter conditions (RCP4.5 and RCP8.5; 1951–2100). The results show that annual precipitation increases of >+150 mm/yr or winter precipitation increases of >+25 mm will lead to significantly higher GWR. GWR is expected to decrease if the precipitation changes are lower than these thresholds. Significant GWR changes are produced only when the temperature change exceeds +2 ◦C. Temperature changes of >+4.5 ◦C limit the GWR increase to +30 mm/yr. This work provides useful insights into the regional assessment of future GWR in cold and humid climates, thus helping in planning decisions as climate change unfolds. The results are expected to be comparable to those in other regions with similar climates in post-glacial geological environments and future climate change conditions.

**Keywords:** groundwater recharge; climate change; thresholds; seasonality; spatiotemporal variations; regional-scale; long-term; HydroBudget model; cold and humid climates; Quebec (Canada)

#### **1. Introduction**

Climate change is already impacting the water cycle everywhere around the world because of precipitation changes and warming temperatures [1]. In particular, it is impacting surface water and groundwater systems in cold and humid climates due to high rates of precipitation change and warming temperatures [2–5]. Because changes at the surface propagate to aquifers through groundwater recharge (GWR), they could have major impacts on groundwater use for human consumption, industry, and agriculture, as well as on groundwater-dependent ecosystems [6–10]. Although the impacts of climate change on groundwater are increasingly studied, the uncertainty associated with simulations of future climatic conditions remains high [9,11–14]. This is even more remarkable in cold and humid climates, where precipitation changes are uncertain (increase or decrease) and where seasonal snow coverage, which leads the annual hydrological cycle, is particularly sensitive to cold season temperatures [2]. A literature review of climate change impacts on

**Citation:** Dubois, E.; Larocque, M.; Gagné, S.; Braun, M. Climate Change Impacts on Groundwater Recharge in Cold and Humid Climates: Controlling Processes and Thresholds. *Climate* **2022**, *10*, 6. https://doi.org/10.3390/cli10010006

Academic Editors: Mohammad Valipour, Sayed M. Bateni and Ying Ouyang

Received: 15 November 2021 Accepted: 4 January 2022 Published: 12 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

groundwater systems in eastern Canada highlighted the wide variability of simulation results from 22 studies spanning the Canadian provinces of Ontario, Quebec, New Brunswick, Nova Scotia, and Prince Edward Island [15]. Using different hydrological, hydrogeological, or integrated models, different downscaling techniques, or different time horizons thus adds further uncertainty to the analysis [15–17], making it difficult to compare future outlooks. Nevertheless, simulation of future conditions remains an essential tool for long-term groundwater resource management in a climate change context.

In cold and humid regions, the geomorphology has largely been shaped by the latest glaciation cycle, groundwater levels are often close to the surface, and unconfined aquifers are generally fed through GWR and connected to superficial water bodies [18–21]. Groundwater recharge is constrained during winter by freezing temperatures that reduce the available liquid water and during the growing season by intensive evapotranspiration rates [4,19,22]. Overall, hydrological systems are highly responsive to changes in the water cycle (e.g., spatio-temporal distribution of precipitation, rainfall intensity, snow accumulation, timing of snowmelt), thus propagating climate changes to the regional hydrology [23]. The impact of climate change in cold and humid regions characterized by an important rise of temperature in the future (especially during the winter season) and by uncertain future precipitation conditions [5,24–26], has been widely studied [2,4,27–33]. The decrease in snow water storage is recognized as a leading cause of low summer stream flows [29,34–36]. As winter temperatures increase, snow cover decreases and winter GWR events become more frequent and are associated with increased winter flow rates in rivers, as evapotranspiration remains very low [4,22,30,37]. Spring peaks of flow rates or GWR become subdued as snowdominated hydraulic regimes transition to rain-dominated regimes [4,11,15,29,33,38–41]. Although not yet well understood, these changes can be important for all groundwater users (humans, industries, ecosystems) and thus need to be studied and forecasted to support future water use policies.

Most climate change impact studies identify ranges of changes in hydrologic variables associated with the climate forcing of the climate scenarios used for simulation [13,27,30,39,41,42]. Reineke et al. [17] observed statistically significant changes in global-scale GWR for different warming levels (+1 to +3 ◦C). Similarly, Döll et al. [7] presented a global analysis of additional hydrologic hazard occurrence resulting from +1.5 and +2 ◦C warming using hydrological indicators, including GWR. However, a range of changes can be insufficient to properly adapt infrastructures and policies to future conditions, as climate change signals usually overlap between climate models and RCPs. To overcome this, Crosbie et al. [13] provided data for water management scenarios using a scale of probability that simulated how future GWR would change from the simulated historical GWR at the Australian scale. Kløve et al. [9] suggested the use of indicators to communicate climate simulation results and representative parameters for use in water resource planning. These indicators of future conditions can be derived from winter low flows [27], GWR volumes [43], or water table depths [44], and can be expressed, for example, as functions of Quaternary deposits [27]. Meanwhile, identification of the evolution of precipitation and temperature changes that would lead to noticeable and statistically significant changes in GWR over time has not yet been undertaken. This assessment of the sensitivity of GWR to changes in climate variables, without a specific distinction between different climate forcing scenarios, would facilitate inter-study comparisons and provide simple and accessible indicators of future conditions for water managers [11].

This work aims to provide new insights into these questions. The objectives are (1) to simulate the effect of climate change on potential GWR in cold and humid climates and post-glacial geological environments and (2) to identify controlling processes and thresholds of GWR changes. As a regional-scale case study, this work focuses on future GWR conditions for the southern region of the Province of Quebec, Canada (36,000 km2), where the hydrology is dominated by long and cold winters. A spatially distributed water budget GWR model calibrated over the 1961–2017 period [22] was used with a set of 12 climate scenarios downscaled from global climate models (GCMs) using RCP4.5 (low emissions) and RCP8.5 (high emissions). The results were used to identify the controlling

processes of GWR changes, as well as temperature and precipitation thresholds that lead to significant long-term changes emerging from future climate uncertainty.

#### **2. Data and Methods**

#### *2.1. Study Area*

A detailed description of the study area can be found in Dubois et al. [22] and is summarized hereafter. The study area is located in the province of Quebec (humid continental climate), between the St. Lawrence River and the Canada–USA border, and between the Quebec–Ontario border and Quebec City (35,800 km2) (Figure 1). It is comprised of the watersheds of eight main tributaries of the St. Lawrence River (numbered W1 to W8 from west to east) (Table 1). Watersheds W1, W2, and W4 have 42%, 83%, and 15% of their total areas, respectively, located in the USA. The topography is flat, with low-elevation areas close to the St. Lawrence River and higher elevations in the Appalachian Mountains. The land cover includes agriculture (42%), forest (45%), wetlands (6%), urban uses (5%), and surface water (2%). The annual average temperature varies between 6.5 (W1, west) and 3.9 ◦C (W8, east), with the west–east cooling gradient also being notable during the cold months (December to March, T < 0 ◦C). The annual precipitation ranges between 952 (W1) and 1123 mm/yr (W4), corresponding to an average of 231 (W3) to 142 mm (W7) of vertical inflow (VI; available liquid water = sum of rainfall and snowmelt) during the cold months.

**Figure 1.** Location of the study area and the watersheds within.


**Table 1.** Watershed characteristics, climate, and potential groundwater recharge (GWR) for the 1961–2017 period (from 22).

\* Part of the watershed is located in the USA—the presented values are only for the Quebec part. Cold months = DJFM; Win. = DJF; Spr. = MAM; Sum. = JJA; Fall = SON; VI = vertical inflow (available liquid water, the sum of rainfall and snowmelt).

The study area includes two main geological units, the sedimentary basin of the St. Lawrence Platform and the metasedimentary Appalachian Mountains. The bedrock is unevenly covered with unconsolidated Quaternary sediments from the last glaciation– deglaciation cycle and is mainly composed of thin and coarse superficial materials deposited on the bedrock in the uphill areas, thick and mixed-grain size deposits in the valleys, and clay covering sandy materials close to the St. Lawrence River. Regional fractured bedrock aquifers flow from the Appalachians to the St. Lawrence River (south/southeast to north/northwest). The aquifers are moderately productive and are in unconfined conditions upstream and semi-confined to confined conditions in the valleys and in the St. Lawrence Lowlands [21]. Dubois et al. [22] estimated the average 1961–2017 regional potential GWR to be 139 mm/yr. They identified preferential recharge zones in the Appalachians, in forested areas, and over coarse deposits and outcropping bedrock. The potential GWR increases from west to east (Table 1). Warmer temperatures in the western watersheds (W1 to W3) are responsible for higher winter GWR rates (more VI) and lower summer and fall GWR rates (more actual evapotranspiration, AET) than in the eastern watersheds. The peak of the spring GWR, which is linked to snowmelt in April, dominates GWR in all the watersheds and corresponds to 44% of the annual GWR rates.

#### *2.2. Simulating Groundwater Recharge with Hydrobudget*

The HydroBudget model (HB) is a water budget GWR model developed to compute spatially distributed potential GWR on grid cells of regional-scale watersheds over long periods of time [45,46]. Using the spatially distributed daily temperature, daily total precipitation, and runoff curve number (RCN—a combination of pedology, land cover, and slopes), HB is driven by eight parameters that require calibration (Table 2) and aggregates the output at a monthly time-step (although daily input data are required). For each daily time-step, HB computes the snow accumulation and melt (two parameters, *TM* and *CM*), tests if the soil is frozen (two parameters, *TTF* and *FT*), and partitions the superficial runoff (based on the RCN and with two parameters, *tAPI* and *frunoff*) from the water that infiltrates into a conceptual soil reservoir (one parameter, *swm*). The AET corresponds to the minimum between the potential evapotranspiration calculated using the formula of Oudin et al. [47] and the available water in the soil reservoir. Part of the residual soil reservoir water is mobilized as potential GWR (one parameter, *finf*). Since HB does not simulate percolation out of the unsaturated zone, the potential GWR represents a maximum of GWR that could reach the saturated zone. For simplification, the simulated potential GWR will be hereon referred to as GWR.


**Table 2.** Description of the HydroBudget parameters and calibrated values from Dubois et al. [22].

Assuming that surface watersheds match hydrogeological watersheds and that rivers drain unconfined aquifers, Dubois et al. [22] calibrated HB by comparing the sum of the simulated superficial runoff and simulated GWR to the measured river flow and the simulated GWR to baseflow estimates (regressive filter on river flow time series). A multiobjective automatic calibration procedure was used with time series of river flow rates from 51 gauging stations over the 1961–2017 period and showed that the simulated variables had a small uncertainty (≤10 mm/yr). Therefore, this regionally calibrated model was considered suitable to be used to simulate future GWR over the study area.

#### *2.3. Climate Scenarios*

A subset of 12 climate scenarios was derived from an ensemble of 54 climate simulations provided by 29 global climate models (GCMs) from the Coupled Model Intercomparison Project—Phase 5 (CMIP5) driven by RCP4.5 and RCP8.5 future greenhouse gas concentrations. The 12-member ensemble (Table 3) was created using the k-means clustering method proposed by Casajus et al. [48]. The climate scenario clustering process was constrained by ten criteria: the change in annual mean temperature and precipitation between the 1981–2010 period and the 2041–2070 period (two variables) and the changes in seasonal mean temperature and precipitation between the same periods (eight variables). Nine out of the 12 clusters comprised climate scenarios based on both RCPs. The algorithm selects the climate scenario closest to the cluster centroid as the candidate (not considering their RCP), as it best represents the average condition of future precipitation and temperature of the cluster. However, CanESM2 (CE2) was hand-picked from its respective cluster to include the Canadian GCM. The subset captures approximately 85% of the initial variance of the ensemble of 54 climate simulations. Casajus et al. [48] showed that this method

retains a good representativity of the uncertainty linked to the climate scenarios between an ensemble of 27 climate scenarios and its subset of six.


The 12 selected simulations were bias-corrected to a 1981–2010 reference dataset (Natural Resources Canada gridded observation database) [49,50] and downscaled to the reference 10 km × 10 km resolution using the quantile mapping approach by Mpelasoka and Chiew [51]. With these scenarios, changes in mean annual temperature and annual precipitation between the 1981–2010 and 2041–2070 periods (ΔT and ΔP, respectively) covered most of the combinations of ΔT and ΔP found in the ensemble of 54 climate scenarios (Figure 2a). ΔT ranged between +0.9 (INM, RCP4.5) and +5.0 ◦C (MIC, RCP8.5) and ΔP ranged between +5 (B1M, RCP4.5) and +200 mm (A13, RCP8.5). The change in mean temperature during the cold months (December to November; ΔTCM) was between +1.1 (INM, RCP4.5) and +6.0 ◦C (MIC, RCP8.5) (Figure 2b). The change in precipitation during the cold months (ΔPCM) was between +17 (MIE, RCP4.5) and +100 mm (GF3, RCP4.5). The warming temperature during the cold months led to ΔVICM between +33 (INM, RCP4.5) and +215 mm (MIC, RCP8.5) (Figure 2c).

#### *2.4. Period Comparisons and Significant Changes*

The simulation period was divided into four 30-year periods: 1981–2010, the reference period, also used as the baseline for the bias correction of the climate scenarios, and three future periods, 2011–2040, 2041–2070, and 2071–2100. The same periods were used to present the results of the 96 GWR scenarios (12 scenarios for 8 watersheds). Each 30-year period was compared to the previous one and to the 1981–2010 reference period to observe the simulated range in future GWR and to identify future GWR changes.

**Figure 2.** Changes (Δ) between the reference period (1981–2010) and the 2041–2070 horizon for the 12 selected climate scenarios in the study area in (**a**) annual precipitation (ΔP) as a function of mean annual temperature (ΔT) within the ensemble of 54 climate scenarios (27 RCP4.5 and 27 RCP8.5), (**b**) cold month (December to March) precipitation (ΔPCM) as a function of mean cold month temperature (ΔTCM), and (**c**) vertical inflow during the cold months (sum of liquid precipitation and snowmelt; ΔVICM) as a function of mean cold month temperature (ΔTCM).

Changes in precipitation, temperature, and GWR were determined to be statistically significant based on the Tukey test (*p* < 0.05), comparing the results of each 30-year period and the previous one or between the future periods and the reference period. The sample size in each group varied between 30 values, when monthly or annual variables for each watershed and each scenario were compared, and 360 values for the monthly or annual variables for each watershed (or grid cells) when all the scenarios were compared.

#### **3. Results**

#### *3.1. Climate Changes for the 1981–2100 Period*

The average evolution of annual precipitation and mean temperature from the 12 scenarios for each watershed and for each 30-year period shows a constant increase throughout the century (Figure 3). All increases between each 30-year period and the previous one (30 years and 12 scenarios corresponding to 360 values for each period) were statistically significant. The range of precipitation and temperature changes, represented by the difference between the minimum and the maximum of the 12 scenarios for each year, increased remarkably from the 1981–2010 period to the 2071–2100 period (Figure 3).

**Figure 3.** Thirty-year period change of mean annual precipitation (**a**) and annual temperature (**b**) of the 12 climate scenarios for the eight watersheds (W1 to W8).

#### *3.2. Groundwater Recharge for the 1981–2100 Period*

The previously calibrated HB model [22] was used to simulate future GWR under the 12 climate scenarios for the entire 1951–2100 period with a monthly time-step and a 500 m × 500 m spatial resolution. Although simulations were performed for the 1951–2100 period, GWR changes were only compared between three future periods (2011–2040, 2041–2070, 2071–2100) and the previous period (including the 1981–2010 reference period). The 30-year moving averages for all GWR scenarios ranged between 126 (W1) and 183 mm/yr (W6) (Figure 4). The ensemble mean GWR increased between +5 (W8) and +17 mm/yr (W1) from the 1981–2010 period to the 2071–2100 period, with maximum increases of +5 mm/yr between two consecutive 30-year periods (Table 4). Six climate scenarios produced GWR rates higher than the ensemble mean (A13, BNU, CMS, GIR, MIC, MRE, RCP8.5, and RCP4.5). The other six climate scenarios (A10, B1M, CE2, GF3, INM, MIE, and RCP4.5) produced GWR rates lower than the ensemble mean. Climate scenarios based on RCP8.5 produced higher GWR rates (although not always the highest).

**Figure 4.** Thirty-year moving average of groundwater recharge (GWR) simulated with the 12 climate scenarios and significant changes (Tukey test, *p* < 0.05) between subsequent 30-year periods (2011– 2040, 2041–2070, and 2071–2100) for (**a**) W1 to (**h**) W8.


**Table 4.** Thirty-year evolution of mean groundwater recharge (mm/yr), range of the ensemble changes (in brackets), and evolution of the cold month groundwater recharge from December to March (T < 0 ◦C) and that from May to November (T > 0 ◦C) (mm) of the 12 climate scenarios for the eight watersheds (W1 to W8).

CM = sum of groundwater recharge for the "cold months", from December to March (T < 0 ◦C); WM = sum of groundwater recharge for the "warm months", from May to November (T > 0 ◦C).

The range of changes in the GWR scenarios was smallest for the 1981–2010 period (2011–2040 for W1 and W2), with values between 22 (W1 and W2) and 28 mm/yr (W4, W6, and W7) (Table 4). It increased markedly in the 2041–2070 period, with values between 42 (W1 and W2) and 60 mm/yr (W5, W6, and W7). It increased even more during the 2071–2100 period, reaching values of between 60 (W1 and W2) and 78 mm/yr (W5). This larger range of the results was due to the increasing range in precipitation changes between the scenarios in the second half of the 21st century (Figure 3).

The climate changes associated with each significant ΔGWR between 30-year periods (not using the 30-year moving average) are reported in Table 5. Although there was a general increase in temperature between 1981 and 2100, a relatively small number of significant inter-period ΔGWR were observed (Figure 4). This could be due to the combined effect of increased evapotranspiration triggered by higher temperature and increased precipitation. As such, the direction of the ΔGWR change was not directly linked to the change in precipitation (ΔP). For example, ΔGWR > 0 (increase) was associated with ΔP < 0 (CMS for the 2071–2100 period, compared to 2041–2070 for W3 and W8), while ΔGWR < 0 (decrease) was simulated with ΔP > 0 (CE2 and GF3 for the 2071–2100 period compared to 2041–2070 for W3 and W8; MIE for the 2041–2070 period compared to 2011–2040 for W3, W6, and W8). An average temperature change (ΔT) of +1.2 ◦C was associated with ΔGWR < 0 (between +0.7 and +2.3 ◦C), while an average ΔT of +1.8 ◦C was associated with ΔGWR > 0 (between +0.2 and +2.8 ◦C). The four climate scenarios based on RCP8.5, representing mainly very humid future conditions, produced statistically significant ΔGWR > 0 for one of the last two future periods, except for A13 in the eastern watersheds (W5 to W8; Figure 4e–h). The climate scenarios based on RCP4.5, representing both moderately and very humid future conditions, produced both ΔGWR < 0 and ΔGWR > 0 for different periods. In addition, the changes between the 1981–2010 and the 2011–2040 periods were not statistically significant (only one significant change in W1 for one scenario; Figure 4a).

Overall, the GWR simulations showed minor variation prior to 2041, and the main changes occurred during the last two future periods, 2041–2070 and 2071–2100. Therefore, only these two future periods will be considered in the rest of the analysis.


**Table 5.** Mean annual precipitation (ΔP) and temperature (ΔT) changes between 30-year periods associated with significant changes in groundwater recharge (ΔGWR) for the 12 climate scenarios and eight watersheds (W1 to W8) (cell color represents the direction of the ΔGWR: orange for decrease and blue for increase, cells are empty when ΔGWR was not significant).

\* Represents RCP8.5 (scenarios without an \* represent RCP4.5). \*\* Period 1 is 2011–2040, Period 2 is 2041–2070, and Period 3 is 2071–2100.

#### *3.3. Inter-Annual Changes in Groundwater Recharge*

The ΔGWR between the two future periods and the reference period for each climate scenario can be represented as a function of different variables (Figure 5). For each scenario, significant ΔGWR < 0 was associated with ΔP < +150 mm for all watersheds except W1 (Figure 5a,b). Inversely, significant ΔGWR > 0 was obtained when ΔP > +150 mm. The few scenarios with ΔP < 0 mm always led to ΔGWR < 0 (some significant, some not; Figure 5b). All the significant ΔGWR < 0 were simulated for +3 ◦C < ΔT < +5 ◦C, while significant ΔGWR > 0 were obtained for +2 ◦C < ΔT < +8 ◦C (Figure 5a,c). ΔGWR seemed to plateau at approximately +30 mm for both ΔT > +4.5 ◦C (Figure 5b,c,h, note triangle markers) and ΔTCM > +6 ◦C (Figure 5f, December to March). Using ΔTCM and ΔPCM showed that ΔGWR < 0 occurred with ΔPCM < +25 mm and +3 ◦C < ΔTCM < +5 ◦C, except for one scenario in W3 and W8 (Figure 5d–f). All significant ΔGWR were simulated with ΔTCM > +3 ◦C (Figure 5f). Significant ΔGWR < 0 were systematically associated with change in cold month GWR (ΔGWRCM) < +25 mm and inversely for ΔGWR > 0 (Figure 5g). Significant ΔGWR < 0 were simulated with scenarios of limited changes in annual simulated AET (+50 mm < ΔAET < +95 mm), while ΔAET > +120 mm were associated with limited GWR increase (plateau around +30 mm) and a ΔT > +4.5 ◦C (Figure 5h). All significant annual ΔGWR were < −15 or > +15 mm (Figure 5b,c,e–h).

Of the 96 GWR simulations, 20 produced a statistically significant ΔGWR between the 1981–2010 and 2041–2070 periods, including 11 based on RCP8.5, and 39 between the 1981–2010 and 2071–2100 periods, including 21 based on RCP8.5 (Table 6). Although scenarios based on RCP8.5 did not always produce significant ΔGWR, they were more likely to produce significant ΔGWR than those based on RCP4.5. The greater number of significant changes simulated for the 2071–2100 period in comparison to the 2041–2070 period confirmed that GWR was more affected with more pronounced climate changes, be it through greater emissions or longer progression.

**Figure 5.** Changes in annual groundwater recharge (ΔGWR) between the reference period (1981–2010) and the 2041–2070 and 2071–2100 periods as a function of (**a**) changes in mean annual temperature (ΔT) and annual precipitation (ΔP), (**b**) annual precipitation changes (ΔP), (**c**) mean annual temperature changes (ΔT), (**d**) mean cold month temperature changes (ΔTCM) and cold month precipitation changes (ΔPCM), (**e**) cold month precipitation changes (ΔPCM), (**f**) mean cold month temperature changes (ΔTCM), (**g**) cold month groundwater recharge changes (ΔGWRCM), and (**h**) annual actual evapotranspiration changes (ΔAET).

**Table 6.** Number of simulations with significant changes in groundwater recharge between the 1981–2010 reference period and the future periods for the eight watersheds (W1 to W8); the number of scenarios based on RCP8.5 producing significant changes is indicated in brackets.


#### *3.4. Spatial Changes in Groundwater Recharge over Time*

The changes in future GWR (future periods vs. the reference period) were analyzed spatially on a cell-by-cell basis with the ensemble of scenarios (Figures 6 and 7). For the months of January, February, and March, and to a lesser extent for the month of December, significant ΔGWR > 0 was simulated for all watersheds, between +1 and >+5 mm for the 2041–2070 period, and mainly >+5 mm for the 2071–2100 period, as well as between +1 and +5 mm for December for the two periods. Although half of the changes were not significant in April for the two future periods, a clear pattern appeared during that month, with −5 mm < ΔGWR < −1 mm in the western portion of the study area and +1 mm < ΔGWR < +5 mm in the eastern portion. In May and June, significant ΔGWR < 0 was simulated, which was lower eastward and for the 2071–2100 period (locally < −5 mm). Generalized significant decreases of −5 mm < ΔGWR < −1 mm were simulated for July, August, and September for the two future periods. The ΔGWR was mainly between −5 and −1 mm from July to November for the two future periods. Non-significant ΔGWR < 0 was simulated in these months in the western and central portions of the study area. Significant ΔGWR < −5 mm was also simulated in the eastern portion for the two future periods in October and to a lesser extent in November.

**Figure 6.** Spatial changes in average monthly groundwater recharge (GWR) between the reference period (1981–2010) and the 2041–2070 period for the 12 climate scenarios.

**Figure 7.** Spatial changes in average monthly groundwater recharge (GWR) between the reference period (1981–2010) and the 2071–2100 period for the 12 climate scenarios.

#### *3.5. Monthly Groundwater Recharge Changes over Time*

The watershed-scale monthly GWR for each period showed significant ΔGWR > 0 simulated for the eight watersheds in December, January, February, and March, with significant increases from 2041–2070 to 2071–2100 from January to March (Figure 8). The range of the ensemble changes for these months also increased remarkably in the future periods in comparison to the reference period. In April, the GWR changes were smaller and the range of the ensemble was smaller. They were mainly significant in the watersheds that are partially located in the USA (W2 and W4). The future GWR in January, February, and March exceeded that of April, which exhibited a peak during the reference period. This can already be noted in the western watersheds (W1 to W4) for the 2041–2070 period and reached similar values in the eastern watersheds (W5 to W8) in March of the 2071–2100 period. Significant ΔGWR < 0 were simulated in May and June for the two future periods and between the two future periods for all watersheds except W1 and W2. Significant ΔGWR < 0 were also simulated in July, August, September, and October for the eight watersheds and between the two future periods in October for the western watersheds (W5 to W8). From May to October, the range of changes of the ensemble was clearly smaller for all watersheds when comparing the reference period with the future periods. While the future GWR was close to zero as early as June and as late as October for the two most western watersheds (W1 and W2), the future GWR reached near-zero values between July and September in the other watersheds. Finally, significant ΔGWR < 0 was simulated for

the eight watersheds in November, again with a smaller range of changes than during the reference period.

**Figure 8.** Monthly groundwater recharge (GWR) for the reference period (1981–2010) and the two future periods (2041–2070 and 2071–2100) for (**a**) W1 to (**h**) W8.

The sum of GWR from December to March increased by +32 mm on average (mean of the ensemble of scenarios), from +27 mm in W1 to +36 mm in W7 and between the 1981–2010 and 2071–2100 periods (Table 4). The sum of the GWR from May to November (months withT>0 ◦C) decreased by −19 mm, from −9 mm in W1 to −26 mm in W8 and between the 1981–2010 and 2071–2100 periods. These seasonal changes were within the range of uncertainty of the annual GWR.

#### **4. Discussion**

#### *4.1. Future Groundwater Recharge Dynamics*

For all watersheds except W1, the GWR for the 2011–2040 period was not statistically different from that of the 1981–2010 period. Of the eight watersheds, significant GWR changes occurred with two to four climate scenarios between the 2011–2040 and 2041–2070 horizons and with four to seven climate scenarios between the 2041–2070 and 2071–2100 horizons (Figure 4). For this reason, the results were compared only between the 2041–2070 and 2071–2100 periods and the 1981–2010 period.

The simulations showed both increases and decreases in GWR in the future, hence markedly increasing the range of possible future conditions from those simulated in the reference period (Figure 4). The climate scenarios based on RCP8.5 were the wettest (140 mm < ΔP < 220 mm, Figure 2) and thus produced increasing GWR rates in the future. Other studies have observed a wide range of hydrological responses to climate change in cold and humid regions or regions with snow-dependent hydrology [2,11,14]. Kurylyk and MacQuarrie [38] simulated increased future annual GWR under four climate scenarios and decreased GWR under three climate scenarios in New Brunswick (eastern Canada). Guay et al. [41] have shown, based on the simulation of future river flows in 305 watersheds in Quebec under 87 climate scenarios, that it was unclear whether future annual river flows would increase or decrease by the 2041–2070 period. Inversely, Sulis et al. [52] mainly simulated decreased future GWR in part of the Chateaugay River watershed (W1) with the integrated CATHY model for the 2041–2065 period (increase under one scenario). These authors used 12 climate scenarios based on the high-emission SRES A2 greenhouse gas projections, with annual precipitation increases of close to 0 to +20% between the future and the reference periods. Differences in future GWR depend on the choice of future horizon and emission scenarios, as well as on the type of model used to derive the GWR [15].

The analysis of monthly recharge allowed major shifts in the intra-annual changes in future GWR to be identified. The results show that winter GWR could significantly increase due to warmer winters and lead to an earlier spring GWR peak. Other studies in eastern Canada have obtained similar results [30,33,37–39]. Similarly, Grinevskiy et al. [53] simulated GWR with an unsaturated zone model (HYDRUS-1D) in 22 sites spread over western Russia (humid climate, cold in the north, temperate in the south) and observed increased GWR during winter, which was linked to wetter and warmer winters in the North, but not in the South of the study area. Such results are reported for cold and humid climates and regions of snow-dominated hydrology with more available liquid water during winter, which is linked to warmer temperatures that affect not only GWR, but the entire hydrologic dynamic [2,11,27,29,34,35,41,54]. These future GWR conditions are supported by observations of past groundwater level time series showing a similar shift in the GWR peak from spring (snowmelt) to winter (rain) in Fennoscandia (Northern Europe, transition between temperate and cold climates) associated with a warming climate between the 1980–1989 and 2001–2010 periods [4].

In the present study, the GWR scenarios showed a statistically significant decrease from May to November (Figure 8). The future GWR was close to zero from July to September (similarly to the reference period), except in the western and warmest watersheds, W1 and W2, where the low flow period began a month earlier (June) and ended a month later (October). Similar results were obtained in different cold and humid climates or regions of snow-dominated hydrology. Guay et al. [41] noticed small or negligible changes in river flows during the summer, a period of the year where flow rates are already very low, extending until October. Expected dryer summer low flow rates were also reported by Addor et al. [11] and Arnoux et al. [27] for the Swiss Alps and by Dieraurer et al. [34] for

watersheds across the Rockies (western North America), and they were linked to reduced snowpack, leading to limited snowmelt contribution to spring flows. In this study, the GWR scenarios were similar for all watersheds during summer, thus increasing the certainty of the expected summer decrease. Despite the high uncertainty in simulations of future hydrologic conditions, Addor et al. [11] reported a good convergence of results toward lower summer flow (90% of the scenarios) for Swiss alpine watersheds, similar to that of the current study. However, Arnoux et al. [27] showed that post-glacial Quaternary deposits can contribute to the mitigation of the impact of climate change on summer low flows in alpine catchments due to their water storage capacity, which supports river low flows during long dry spells. This is an indication that water-bearing unconsolidated superficial materials could be an indicator of watershed response to climate change.

Aygün et al. [2] showed that the hydrology of cold and humid regions (northern regions of North America and Eurasia outside of the permafrost zone) with near-freezing annual temperatures were more sensitive to climate change than regions with substantially colder climates. The current study showed differences in the watershed response from west to east that followed the regional temperature gradient (decrease of mean annual temperature). These findings were most likely possible because of the use of a single model across the region and a robust knowledge of the past dynamics. Larocque et al. [15] did not find such a clear trend from west to east in their review of modeling studies of climate change impacts on groundwater systems in eastern Canada.

#### *4.2. Climate Changes Impacting Groundwater Recharge*

The groundwater recharge changes became statistically significant when ΔGWR was < −15 or > +15 mm for the two future periods (Figure 5). More specifically, small GWR changes could not be interpreted as being different from the simulated variability of the 1981–2010 reference period for one to five of the 12 scenarios for the 2041–2070 period and three to six of the scenarios for the 2071–2100 period (Table 6). The increasing number of scenarios with significant changes for the 2071–2100 period is coherent with results from Goderniaux et al. [42] in the Geer Basin (Belgium), where projections of groundwater levels obtained using an ensemble of 30 climate scenarios became greater than the variability of the 1961–1990 period only in 2085. Similarly, using an ensemble of 54 climate scenarios, Addor et al. [11] demonstrated that flow rate changes in alpine catchments became significantly different from those of the 1980–2009 reference period only after the 2050 horizon. They showed that significant changes were simulated even under the climate scenarios with the lowest emissions based on RCP2.6 (not used in this study). In contrast, this study showed that climate scenarios based on RCP8.5 did not systematically produce significant changes between the future periods and the reference period, although they tended to simulate significant changes and higher future GWR than scenarios based on RCP4.5 more often. Henceforth, using a large ensemble of climate scenarios appears to be necessary to provide a representative sample of possible future precipitation and temperature.

One of the main novelties of this work lies in the identification of climate conditions leading to statistically significant changes in future GWR. Significant ΔGWR < 0 was simulated only with ΔP < +150 and ΔPCM < +25 mm. ΔP < 0 always led to ΔGWR < 0, but the latter was not necessarily significant. Inversely, ΔGWR > 0 were significant only with ΔP > +150 and ΔPCM > +25 mm. Therefore, ΔP ≈ +150 and ΔPCM ≈ +25 mm appear to be regional thresholds for determining the direction of future GWR changes.

Another contribution of this work was to determine that significant ΔGWR < 0 was systematically associated with ΔT and ΔTCM ranging between +3 and +5 ◦C, while significant ΔGWR > 0 was found for +2 ◦C < ΔT < +8 ◦C and +3 ◦C < ΔTCM < +11 ◦C. Interestingly, ΔT > +4.5 ◦C (or ΔTCM > +6 ◦C) led to ΔAET > +120 mm, thus limiting ΔGWR to +30 mm. Therefore, ΔT ≈ +2 ◦C and ΔTCM ≈ +3 ◦C appear to be regional thresholds for significant GWR changes (increase or decrease), while ΔT > +4.5 ◦C triggers GWR increase. These temperature thresholds control future GWR through the modification of the cold month

hydrology and the evolution of ET, which also depends on the adaptation of vegetation to climate change.

This study demonstrated that, on an annual basis, ΔGWRCM > +25 mm compensated for decreased GWR during the rest of the year and produced statistically significant ΔGWR > 0. This is coherent with numerous previous studies showing that the seasonality of the entire hydrologic dynamic (GWR, groundwater storage, groundwater level, stream flow) in cold and humid climates or regions of snow-dependent hydrology was affected by the increase in available liquid water during warmer winters, counterbalancing the decreasing availability of water during summer [2,4,11,27,30,35,39,41]. Rivard et al. [39] observed that changes in future GWR were most sensitive to winter temperature in simulations with a spatialized water budget model in Nova Scotia (HELP, eastern Canada, not overlapping the current study area). This was more due to increased amounts of liquid winter precipitation that was readily available for infiltration than to changes in precipitation amounts. Interestingly, Wright and Novakowski [33] showed that winter recharge events on a fractured bedrock (Ontario, Canada, frozen during winter) could bypass the frozen soil and reach unfrozen fractures at the soil/bedrock interface. They concluded that winter rainfall events could produce more GWR than during the rest of the year, thus making the precipitation form and amount during this period a sensitive GWR variable. The differences in these studies may be due to their respective scales. The local scale associated with GWR estimates based on well observations used in some studies [4,33] is in dire contrast to the 250 m × 250 m resolution used by Rivard et al. [39] or the 500 m × 500 m in the HydroBudget model. In addition, the sub-hourly sampling time-step of other studies [33,55] is not comparable to the daily time-steps aggregated into monthly inter-annual results presented here. Nevertheless, all the available studies for Eastern Canada confirm the importance of future winter GWR in the overall annual GWR dynamic, as well as the importance of capturing local-scale (meter order) processes in regional-scale GWR simulations.

From a different perspective, Sulis et al. [52] showed that changes in future GWR in a sub-watershed of W1 were linked to intra-annual patterns of the climate scenarios (more snowmelt during winter, less rain during the fall, the duration of successive days with daily precipitation > 1 mm/d) rather than being related to annual precipitation changes. The integrated CATHY model (daily time-step) seemed sensitive to the dryness conditions of the soil [56], thus inducing more percolation through the unsaturated zone (GWR) for climate scenarios with regular summer rainfall events than for scenarios with more intense but less frequent rainfall events. Similar conclusions were reached by Wright and Novakowski [33] at the well scale in a fractured bedrock aquifer for winter GWR events in Ontario. Finally, Rathay et al. [55] observed that increasing rainfall intensity, from <1 mm/h to > 1 mm/h, produced a decrease in the rainfall–groundwater level cross-correlation coefficients in a bedrock aquifer in the temperate climate of British Columbia (Canada). Although they did not identify a rainfall intensity threshold limiting GWR, these authors concluded that more intense rainfall events produced more surface and subsurface runoff rather than increasing GWR rates. Although these studies highlighted that precipitation intensity can be an important factor for future GWR changes in humid climates, the sensitivity of GWR to this parameter was not a focus of the current study.

#### *4.3. Future Groundwater Recharge Simulation in Cold and Humid Climates*

The clustering method used to select the subset of climate scenarios was based on ten criteria including changes in seasonal and annual precipitation, as well as changes in temperature, but did not include changes in precipitation intensity. Although recent work has projected the intensification of year-round precipitation in North America [57], precipitation intensity changes for the province of Quebec are not yet clear [25]. Further research needs to assess the impact of this variable on increasing or decreasing future GWR in cold and humid climates on intra-annual and inter-annual time scales.

Considering the range of changes in future recharge, understanding GWR under future conditions probably lies mostly in the capacity to adequately simulate GWR during the cold months, the period corresponding to the greatest changes in terms of absolute value in the study area. Stream flow or GWR simulations in cold and humid climates are sensitive to snow-related calibration parameters, such as the melting temperature and melting coefficient [22,58]. However, Melsen and Guse [59] showed that these parameters were less sensitive when simulating river flow in 605 USA watersheds under future conditions with decreased snowpack. Therefore, an evolution of the snow-related parameters could be expected under future conditions. Improving the simulation of winter GWR in cold regions will necessitate a better understanding of the roles of snow dynamics and soil frost in changing conditions, and future work should be aimed at calibrating these parameters for long-term regional-scale simulations.

The current study was based on the HB model, which was calibrated and validated over an exceptionally long period of time (57 years from 1961 to 2017), ensuring satisfying representativeness of the long-term and regional-scale hydrological dynamics [22]. The resulting GWR scenarios used constant model parameters over time under the hypothesis that the system was stationary in time and no significant land-use change occurred. However, Jaramillo et al. [60] linked a 40-year increase in AET rates of 65 mm/yr in the Stockholm region (Sweden, temperate to cold climate transition zone) to land-use change, with the massive conversion of semi-natural grasslands (mowing) to cereal and fodder harvesting at the beginning of the 20th century. For Sweden as well, Destouni et al. [61] compared the evolution of evapotranspiration (ET) and runoff (R) for nine watersheds in temperate and cold climates that remained stationary in time or were affected by hydropower and non-irrigated agricultural development during the 20th century. Despite precipitation and temperature increases, they found that ET and R remained stable in unregulated watersheds, while hydropower development increased ET and decreased R, and agriculture development increased both ET and R. These hydrological changes impact the regional water budget, and therefore most likely propagate to GWR. Alternatively, Guerrero-Morales et al. [62] found that land cover changes accounted for 25% of the GWR decrease in an urbanized watershed in western Mexico (warm and humid climate) under climate change conditions by the 2050 horizon. Although Kløve et al. [9] and Taylor et al. [10] stated that climate change studies should consider land-use change, integrating land-use scenarios into future GWR simulations in cold and humid climates has not been widely reported in the scientific literature. Further study of this important question could lead to the identification of other factors than climate that determine the extent of possible GWR changes. Considering land-use change would also probably increase the uncertainty of future GWR simulations [63].

Reinecke et al. [17] concluded on the importance of coupling biosphere dynamic simulations to long-term GWR simulation, especially at the global scale, where increases in atmospheric CO2 concentrations could lead to more active vegetation, which would, in turn, impact GWR estimates. Koirala et al. [64] showed that vegetation had a large impact on the water budget through AET, especially in humid climates. To avoid using scenarios of future solar radiation and other climate variables that are less readily available, more difficult to bias-correct, and may introduce additional uncertainty compared to the more common temperature and precipitation scenarios [25], the maximum daily ET in HB was based on the simple formula from Oudin et al. [47], which only used daily temperature, latitude, and Julian day (as a proxy for extraterrestrial radiation). However, considering the regional scale of the study and the long-term simulation period, more work should be dedicated to improving AET simulations for cold and humid climates, especially considering the uncertainty related to plant adaptation to warmer climates. This could impact the temperature thresholds identified in this study. Specific AET calibration could be developed using spatialized time-series of the measured AET, or the impact of coupling biosphere dynamics and GWR at the regional scale could be tested.

#### *4.4. Using These Results for Adaptation*

Studies clearly show that GWR in cold and humid climates could follow different paths of change depending on specific climate conditions, geology, morphology, and land use [15,38]. Taking into consideration uncertainty in future climate conditions is another major challenge, as this study showed with future GWR that can either increase or decrease at the regional scale depending on the climate scenario It is thus extremely difficult to provide concrete recommendations to water managers despite the increasing body of knowledge [14].

Nevertheless, several patterns in the future evolution of GWR emerged with a relatively high level of confidence. For example, the significant projected decrease in GWR from May to November as soon as the 2041–2070 period and the substantial increase in GWR from December to March clearly stand out. A cold month GWR increase of > +25 mm will compensate for the decrease throughout the rest of the year, suggesting stable groundwater resources. Additionally, this work provides threshold values for changes in precipitation and temperature that lead to likely increases or decreases in future GWR (Figure 9). These thresholds could be used in integrated water resource management plans, where they could trigger specific actions (e.g., if local warming reaches 1.5, 2, or 3 ◦C, associated with stable precipitation increase of +50 or +100 mm). Although they would probably be similar in other cold and humid climates in post-glacial geological environments, these thresholds will need to be tested in different contexts.

**Figure 9.** Annual groundwater recharge changes (ΔGWR) between the reference period (1981–2010) and future periods (2041–2070 and 2071–2100) for the eight watersheds (W1–W8) and 12 scenarios. The associated precipitation and temperature thresholds are displayed on the right, and the gray zone indicates the −15 to +15 mm non-significant change range in GWR. CM stands for cold months (December to March).

In cold and humid climates, GWR generally represents the actual aquifer renewal rates—the total flow discharging to superficial water bodies [4,19,21]. A decrease in future GWR from May to November means that groundwater inflow into superficial water bodies and groundwater levels will decrease when water demand for drinking water, agriculture, industrial purposes, hydroelectricity, and recreation is the highest [19,65] and when river flows come almost exclusively from a connected aquifer. Considering the high confidence

in the simulation of the decreasing GWR from May to November, it is expected that water use conflicts will increase in future decades.

The identified thresholds were related to potential GWR, i.e., the maximum GWR that can reach the saturated zone [22]. Future changes in actual GWR are expected to be closely linked to those in potential GWR. For scenarios and periods where potential GWR is expected to decrease, actual GWR will most likely decrease as well due to the reduction of available water. Inversely, for expected increases in potential GWR, actual GWR changes would vary depending on the AET rates. Future work studying the propagation of these changes should focus on the periods of expected potential GWR increases.

#### **5. Conclusions**

In cold and humid climates, the impact of climate change will propagate in groundwater systems and more broadly to regional hydrologic dynamics through GWR. Estimates of changes in GWR under future climate conditions are therefore strategic for long-term water resource management. This work has provided new data for assessing climate change impacts on GWR and to identify controlling processes and thresholds for cold and humid climates. One of the outcomes was the simulation of the first set of 12 transient regionalscale GWR scenarios for the 1951–2100 period in southern Quebec. Simulated using a water budget model and a set of 12 climate scenarios maximizing the future climate variability (12 GCMs using RCP4.5 and RCP8.5), the spatio-temporal GWR scenarios showed notable changes occurring in the 2041–2070 and 2071–2100 periods. Warming temperatures were between +1 and +5 ◦C at the 2041–2070 horizon (in comparison with the 1981–2010 reference period), and the precipitation change pattern was more variable, including an increase of +10% to +80% in the available liquid water between March and December. Increasing and decreasing annual GWR was simulated. However, major impacts were found in the monthly dynamics, with a statistically significant decrease in future GWR from May to November compensated by a statistically significant increase in future GWR from December to March. The periods of null or very low GWR rates were lengthened by one month in June and October for the warmer watersheds. Overall, the average annual GWR change was positive if the increase in future cold month GWR was higher than +25 mm, offsetting the decrease for the rest of the year. Such results were coherent with previous findings in other regions of cold and humid climates.

The novelty of this work lies in linking changing climate conditions to the direction and amplitude of statistically significant changes in future regional GWR through specific precipitation and temperature change thresholds. All significant changes in GWR were >+15 or <−15 mm/yr and were only produced by warming temperatures >+2 ◦C. A significant decrease in future GWR was always simulated under future increases in annual precipitation of <+150 mm and cold month precipitation changes of <+25 mm, along with warming temperatures of between +3 and +5 ◦C (for annual and cold months). A significant increase in future GWR was systematically simulated under increases in annual precipitation of >+150 mm and cold month precipitation increases of >+25 mm, along with warming temperatures of >+2 ◦C. A future temperature increase of >+4.5 ◦C produced more intense AET rates, thus limiting the increase in future GWR to approximately +30 mm, irrespective of the precipitation increase. These thresholds are sufficiently straightforward for general use and for integrated water resource management plans.

**Author Contributions:** All authors contributed to writing the manuscript. E.D., M.L., and S.G. developed the approach. M.B. selected climate scenarios through cluster analysis. Simulations and figure preparation were done by E.D., M.L. obtained the research grant and supervised the research. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Quebec Ministry of Environment and Climate Change (Ministère de l'Environnement et de la Lutte contre les changements climatiques).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The simulated scenarios presented in this work are available here: https://doi.org/10.5683/SP3/SWH4O1.

**Acknowledgments:** The authors are grateful to the Ouranos Consortium for providing downscaled climate scenarios and acknowledge the model output data from the World Climate Research Programme's Coupled Modelling Intercomparison Project Phase 5 (CMIP5) as well as the gridded observation data made available by Natural Resources Canada's (NRCan).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Hydro-Meteorological Trends in an Austrian Low-Mountain Catchment**

**Gerald Krebs 1,2,\*, David Camhy <sup>1</sup> and Dirk Muschalla <sup>1</sup>**


**Abstract:** While ongoing climate change is well documented, the impacts exhibit a substantial variability, both in direction and magnitude, visible even at regional and local scales. However, the knowledge of regional impacts is crucial for the design of mitigation and adaptation measures, particularly when changes in the hydrological cycle are concerned. In this paper, we present hydrometeorological trends based on observations from a hydrological research basin in Eastern Austria between 1979 and 2019. The analyzed variables include air temperature, precipitation, and catchment runoff. Additionally, the number of wet days, trends for catchment evapotranspiration, and computed potential evapotranspiration were derived. Long-term trends were computed using a non-parametric Mann–Kendall test. The analysis shows that while mean annual temperatures were decreasing and annual temperature minima remained constant, annual maxima were rising. Long-term trends indicate a shift of precipitation to the summer, with minor variations observed for the remaining seasons and at an annual scale. Observed precipitation intensities mainly increased in spring and summer between 1979 and 2019. Catchment actual evapotranspiration, computed based on catchment precipitation and outflow, showed no significant trend for the observed time period, while potential evapotranspiration rates based on remote sensing data increased between 1981 and 2019.

**Keywords:** hydrological research basin; precipitation; temperature; long-term trends; climate change; evapotranspiration

#### **1. Introduction**

It is well documented that the climate is changing [1–3]. Impacts are seen as globally rising temperatures [2,4] with a reduced number of cold days and nights and an increased number of warm days and nights [4], an altered depth [5–7] and duration of snow and ice cover [6–8], changing precipitation [4,9–11] and river flow regimes [12–14], or an increased number of extreme events [2,4,15]. However, the magnitude and impact direction of major climate variables, such as temperature, precipitation, catchment runoff, and evapotranspiration in both climate observations and projections vary significantly at the global and regional scale [16,17].

While there is a consensus on global warming [2] supported by many studies (e.g., [15,18,19]), some areas experienced decreasing mean, maximum, or minimum temperatures 1951-2002 [20]. Precipitation observations indicate minor global changes despite a large, compensating variability with a decrease observed in the subtropics and Southern Europe [21,22], Southern Asia and Africa and increases observed in North America, South America, Eurasia, and North and Central Europe [11,22–24]. Furthermore, a seasonal shift of precipitation has been reported (e.g., [19,25]).

With respect to catchment runoff, a decrease was observed for some basins in China [18,26], while an increase in runoff was reported for other Chinese basins [27] or North-Eastern USA [28]. Blaschke et al. [24] report only minor changes in runoff for

**Citation:** Krebs, G.; Camhy, D.; Muschalla, D. Hydro-Meteorological Trends in an Austrian Low-Mountain Catchment. *Climate* **2021**, *9*, 122. https://doi.org/10.3390/cli9080122

Academic Editor: Mohammad Valipour

Received: 2 June 2021 Accepted: 26 July 2021 Published: 29 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Austrian catchments in the past 50 years but predict a future runoff reduction for summer and an increase for winter. As for precipitation, a seasonal runoff shift has been reported (e.g., [13]). Furthermore, an increase in flooding events, particularly in Alpine areas, has been observed [14].

Several studies report increasing potential evapotranspiration trends for most of the Northern hemisphere (e.g., [19,29–33]), while China experienced decreasing evapotranspiration rates over the past 50 years [34]. Some of these studies confirm the trend that dry areas become drier and wet areas become wetter, while some contradict [23,35].

The validation of observations is one of the most important tasks during hydrological assessments as faulty data obviously provoke wrong analysis results and conclusions. At the same time, particularly the validation of precipitation measurements is very demanding due to the spatial and temporal variability of rainfall and its stochastic nature. An appropriate validation strategy depends on several factors, such as the spatial distribution of stations, the recording and analysis frequency or the type of measurement device. While there is no standardized procedure that is generally applicable, validation strategies commonly comprise the following steps: (i) identification of documented defects, (ii) device-specific boundaries, (iii) climatological boundaries, (iv) temporal variability, (v) intra-stational validation, and (vi) inter-stational variability [36,37].

The literature shows that the impact of climate change is widely acknowledged. At the same time, it is obvious that the impacts greatly vary at a regional and even local scale. However, this knowledge is crucial to develop measures to mitigate and counteract hydrological climate change impacts at the regional and local scale. Furthermore, we aim to investigate whether large-scale climate observations or projections also hold for smaller catchments where hydro-meteorological conditions may be very site-specific. For this purpose, we analyzed the hydro-meteorological data from a hydrological research catchment in Styria, Eastern Austria, which has been monitored since 1979. Analyzed climate variables include precipitation depth and intensities, number of wet days, air temperature, river flow, and actual and potential evapotranspiration.

#### **2. Materials and Methods**

#### *2.1. Hydrological Research Catchment Pöllau*

The hydrological research basin (HRB) *Pöllau* was established in 1978 [38,39] and is currently operated by the Institute of Urban Water Management and Landscape Water Engineering at Graz University of Technology in cooperation with the Department 14 of the Federal State Styria. The decision to establish an HRB in the *Pöllau* sub-basin was based on a number of reasons: (i) the confining arched mountain ridge allows a clear delineation of the catchment, (ii) the loamy soils are characterized by low storage capacities, minimizing the influence of subsurface flow on catchment hydrology, and (iii) the climate of the catchment with heavy storm events in the summer and relatively dry winters is representative for the Eastern alpine foothills [40]. The catchment covers 58.3 km<sup>2</sup> and is located in Styria, Austria, about 60 km northeast of the city of Graz (Figure 1). The elevation of the catchment ranges from 398 to 1279 m, and the catchment land-cover is dominated by forest (ca. 43.8%) and grass- and cropland (ca. 51.8%) with a low degree of discontinuous urban fabric (ca. 4.4%) [41].

The catchment comprises two main sub-catchments that are monitored: (i) the subcatchment Saifenbach/Dürre Saifen covering 23 km<sup>2</sup> (monitored 1997–2005 and since 2018) and (ii) the sub-catchment Prätisbach covering 21 km<sup>2</sup> (monitored since 1980). Additionally, the discharge at the joint catchment outlet of both the sub-catchments has been monitored since 1980. Characteristic catchment properties are given in Table 1.

**Figure 1.** Overview of the catchment *Pöllau* (discharge measurement *A*) with the sub-catchment *Prätisbach* (discharge measurement *B*) in the west and the sub-catchment *Dürre Saifen* (discharge measurement *D*) in the east, and the locations of the precipitation measurements.



#### *2.2. Data*

The first precipitation measurement gauge in the HRB *Pöllau* was installed in 1979 (*1*, see Figure 1 and Table 2). During the following year (1980), an additional five precipitation gauges were installed, and two stream gauges (the catchment outlet *A* and the sub-catchment *B*) were constructed and taken into operation. The precipitation monitoring at the meteorological station (*7*) started in 1982, whereas the observation of climate variables started in 1991. The stream gauge *C* started operation in 1988 but was destroyed during a massive flood in 1997. The gauge was then reconstructed in 2000, but after further flood damage in 2007, the gauge was not put back into operation. The stream gauge *D* was constructed in 1997, but due to the challenging measurement location, monitoring was abandoned in 2005. The gauge was reconstructed 500 m upstream in 2018 and is, together with the gauges *A* and *B*, currently operating.

The currently operated precipitation gauges are rather symmetrically distributed over the catchment area and located at elevations between 420 and 1040 m.a.s.l. Initially, all 7 precipitation gauges were tipping buckets with a resolution of 0.1 mm. Since the year 2011, 6 stations have been equipped with rain scales (type Ott *Pluvio*<sup>2</sup> [42]) operated at a 1 min recording interval. The rain scales in the catchment are not equipped with heaters, which hampers the record of snowfall. Due to this shortcoming, only snow that stayed on the gauges and thereafter melted was recorded. The currently operated stream gauges monitor the entire catchment outflow (*A*) and the two main sub-catchments (Figure 1). The stream

gauges are equipped with pressure sensors, calibrated with rating curves, and record at a 10–15 min interval.

**Table 2.** Stations, altitude (m.a.s.l), measured variables: *WL* (water level), *WT* (water temperature), *P* (precipitation), *T* (air temperature), *p* (air pressure), *rH* (relative humidity), *Ra* (solar radiation), *ST* (soil temperature), *SM* (soil moisture), *WS* (wind speed), *WD* (wind direction), and data availability.


Land-cover changes due to urbanization, agriculture or forestation can, along with potential climate change, significantly affect the catchment water balance, hampering the attribution of observed long-term changes to a single driver. Therefore the catchment land-use was analyzed based on the CORINE land-cover datasets available for the years 1990, 2000, 2006, 2012, and 2018 [43]. Between 1990 and 2018, forested areas decreased by 0.27 km2 and were mostly replaced by settlements growing by 0.3 km2, while the fraction of agricultural landuse remained almost constant with a decrease of 0.03 km2 (Table 3).

**Table 3.** Land-cover in the catchment *Pöllau* in 1990, 2000, 2006, 2012, and 2018 based on the CORINE land-cover datasets [43].


#### *2.3. Data Validation*

To exclude as much doubtful data as possible from the subsequent analysis, the available measurements were first validated on a daily basis according to the following procedure: (i) identification of documented defects, (ii) device-specific boundaries, (iii) climatological boundaries, (iv) temporal variability, (v) intra-stational validation, and (vi) inter-stational variability [36,37].

The validation steps (i)–(vi) were applied for rainfall and discharge observations. The comparison of daily precipitation observations after validation shows a good correlation (Pearson correlation 0.91), allowing the conclusion that the seven stations mostly recorded similar values (Figure 2). Discharge measurements were validated using cumulative sums of available gauges. An inter-stational validation for temperature data was not directly possible, as this variable is recorded at only one location within the catchment. However, the general observed pattern was compared with regionally available temperature observations for consistency.

**Figure 2.** Scatter of daily recordings of each station against each station (Pearson correlation 0.91).

#### *2.4. Data Analysis*

Long-term hydrological trends and their significance were computed using the nonparametric modified Mann–Kendall test [44] to reduce the influence of serial correlation. Additionally, the Theil–Sen robust estimate was computed [45,46] to evaluate the magnitude of the trend. This approach has been successfully used to assess climate developments in numerous earlier studies (e.g., [47–50]) and was therefore applied in the current study.

Long-term trends of air temperatures were analyzed based on mean annual temperatures, on the one hand, and on mean seasonal temperatures recorded at the climate station *7* on the other hand. The seasons were defined as spring (March, April, May), summer (June, July, August), autumn (September, October, November) and winter (December, January, February). Seasonal trends were computed as annual trends that might be balanced by seasonal changes.

The conducted precipitation analyses comprised long-term trends of annual and seasonal (seasons as defined above), precipitation depths, and long-term trends of precipitation intensities for different durations (60, 120, 240 min). Precipitation depth was analyzed as the catchment mean sum (mean of the station recordings that fulfilled the validation criteria). Additionally, the trend of the number of annual wet days in the catchment was analyzed.

Long-term trends for the catchment discharge were analyzed for gauge *A*, while the remaining gauges were utilized for data validation only.

As for precipitation and temperature, long-term flow trends were also analyzed at a seasonal scale to identify temporal shifts in stream flow behavior.

The catchment water balance was computed based on observed precipitation runoff to assess the long-term development of actual evapotranspiration in the catchment. The computation includes a number of simplifications: (i) groundwater outflow of the catchment is not considered (no data available), (ii) land-cover changes were minor in the catchment during the observation period and therefore not further considered, and (iii) only years are taken into account, where available data allow for computation of annual runoff values. The simplifications yield the following water balance:

$$ET = P - R \tag{1}$$

where *ET* is actual evapotranspiration (mm), *P* is observed catchment precipitation (mm), and *R* is observed catchment runoff (mm). Additionally, potential catchment evapotranspiration (PET) was computed based on remote sensing data from the Copernicus Climate Data Store [51] as global radiation observations at the catchment climate station were insufficiently complete to allow for PET computation. PET was computed according to de Bruin et al. [52] using temperature, surface solar radiation and top of the atmosphere incident solar radiation. This approach was selected as wind speed data recordings were incomplete at the catchment climate station. Wind speed substantially influences PET, and thus, this simplification may underestimate computed PET rates.

#### **3. Results**

#### *3.1. Temperature Trends*

The mean annual air temperature at the climate station *7* between 1991 and 2019 is 9.8 ◦C, with the maximum annual mean recorded in 1995 (11.5 ◦C) and the minimum annual mean recorded in 1991 (7.6 ◦C). The long-term development of the mean annual temperature shows a negative trend with decreasing annual mean air temperature recordings (Figure 3).

While the development of annual minima shows no significant trend, annual temperature maxima were increasing between 1991 and 2019. The mean annual minimum 1991–2019 is −13.4 ◦C, with the lowest recording in 2009 (−19.2 ◦C) and the highest recording in 2015 (−8.3 ◦C). The mean annual maximum 1991–2019 is 32.5 ◦C, with the lowest recording in 1997 (29.1 ◦C) and the highest recording in 2003 and 2016 (37.4 ◦C) (Figure 3).

**Figure 3.** Annual mean (black), maxima (red) and minima (blue) of the air temperature and their trends 1979–2019.

The mean winter temperature between 1991 and 2019 is 0.3 ◦C shows a decreasing trend, with the lowest value recorded in 2009 (−3.2 ◦C) and the highest recording in 1994 (2.7 ◦C). Both winter minima (mean of −13.4 ◦C, with the lowest recording in 2009 (−19.2 ◦C) and the highest recording in 2015 (−8.3 ◦C)) and maxima 1991–2019 (mean of 16.5 ◦C, with the lowest recording in 1996 (10.6 ◦C) and the highest recording in 2011 (20.3 ◦C)) show no significant trend (Figure 4 top left).

The mean spring temperature between 1991 and 2019 is 10.0 ◦C and shows a similarly decreasing trend as observed for winter. The lowest mean was recorded in 2009 (7.6 ◦C) and the highest recording in 1995 (12.1 ◦C). The trend of spring minima is decreasing around a mean minimum of −5.5 ◦C, with the lowest recording in 2018 (−16.3 ◦C) and the highest record observed in 2011 (−1.2 ◦C). The trend of spring maxima 1991–2019 is also decreasing around 27.0 ◦C, with the lowest recording in 1991 (23.1 ◦C) and the highest recording in 1999 (30.9 ◦C) (Figure 4 top right).

The mean temperature during summer between 1991 and 2019 shows no trend staying at 19.9 ◦C, with the lowest recording in 2008 (17.5 ◦C) and the highest recording in 2003 (22.2 ◦C). The trend of summer minima is also not significant at 6.3 ◦C, with the lowest recording in 2006 (2.7 ◦C) and the highest recording in 2019 (12.7 ◦C). The trend of summer maxima 1991–2019 is increasing around 35.5 ◦C, with the lowest recording in 1997 (29.1 ◦C) and the highest recording in 2003 and 2016 (37.4 ◦C) (Figure 4 bottom left).

The trend of the mean autumn temperature 1991–2019 is not significant around 9.8 ◦C, with the lowest recording in 2008 (6.0 ◦C) and the highest recording in 1995 (13.4 ◦C). The trend of autumn minima is decreasing around −4.9 ◦C, with the lowest recording in 2008 (−9.8 ◦C) and the highest recording in 1995 (1.0 ◦C). The trend of autumn maxima 1991–2019 is increasing around 25.7 ◦C, with the lowest recording in 2010 (22.7 ◦C) and the highest recording in 2015 (31.5 ◦C) (Figure 4 bottom right). A comprehensive summary of the observed temperature trends, including statistical trend properties, is given in Table 4.

**Figure 4.** Mean (black), maxima (red) and minima (blue) of the air temperature and the trend 1979–2019 for the winter (**top left**), spring (**top right**), summer (**bottom left**), and autumn (**bottom right**).


**Table 4.** Summary of the climate variable trends for the catchment *Pöllau*.


**Table 4.** *Cont.*

*3.2. Precipitation Trends*

3.2.1. Precipitation Depth

The mean annual precipitation shows no significant trend between 1979 and 2019 around 608.9 mm, with the maximum mean recorded in 2014 (807.2 mm) and the minimum recorded in 2001 (364.3 mm). The annual maximum at a single station was recorded in 1996 at *4* (829.2 mm) and the annual minimum in 2001 at *7* (340.4 mm) (Figure 5). It is to be noted that the circumstance in which the rain scales are not equipped with heaters allows for the recording of only melted snowfall. Thus, winter and partly spring precipitation trends depend both on snowfall and melting occurring. The same is true for derived precipitation intensity trends.

The seasonal precipitation 1979–2019 shows an increasing trend for summer (June, July, August), while no significant trend was detected for spring (March, April, May), autumn (September, October, November), and winter (December, January, February) 1979–2019 (Figure 6).

**Figure 5.** Annual precipitation of the 7 stations (25% and 75% percentile, median (red)) and trend 1979–2019 (dashed line).

**Figure 6.** Precipitation of the 7 stations (25% and 75% percentile, median (red)) and trend 1979–2019 (dashed line) for winter (**top left**), spring (**top right**), summer (**bottom left**), and autumn (**bottom right**).

The mean winter precipitation in the catchment 1979–2019 was 73.3 mm, with the highest recording in 2013 (139.5 mm) and the lowest recording in 1998 (16.3 mm) (Figure 6 top left). The mean precipitation falling in the winter season accounted for 12% of the mean annual precipitation 1979–2019.

The mean spring precipitation accounted for 151.9 mm for 25% of the mean annual precipitation 1979–2019. The largest spring precipitation was recorded in 1985 (272.5 mm) and the smallest in 2003 (66.7 mm) (Figure 6 top right).

The mean summer precipitation shows a clearly increasing trend around 222.0 mm, accounting for 36% of the mean annual precipitation 1979–2019. The largest summer precipitation was recorded in 2018 (416.4 mm), and the smallest value was recorded in 1984 (99.1 mm) (Figure 6 bottom left).

The mean autumn precipitation 1979–2019 was around 168.2 mm, accounting for 27% of the mean annual precipitation. The largest autumn precipitation was recorded in 1993 (273.5 mm) and the smallest precipitation in 2019 (78.7 mm) (Figure 6 bottom right).

#### 3.2.2. Wet Days

The number of wet days in the catchment remained constant, with a mean of 78 wet days per year (Figure 7). The highest number of wet days was recorded in 1979 with 105, while the smallest number of rainfall days was recorded in 2019 with 54.

#### 3.2.3. Precipitation Intensities

Precipitation intensities for a duration of 60 min showed no significant trend at an annual level as well as for summer and autumn. However, an increasing trend was detected for winter and spring 1979–2019. Annual intensities for a duration of 120 min showed no significant trend as well as for winter and autumn, while spring and summer experienced

increasing intensities (Figure 8). The trend for a longer duration of 240 min was not significant for winter and autumn as well as annually. However, as for the duration of 120 min, intensities were increasing for spring and summer. A comprehensive summary of the observed precipitation trends, including statistical trend properties, is given in Table 4.

**Figure 7.** Annual wet days recorded at the 7 stations (25% and 75% percentile, median (red)) and trend 1979–2019 (dashed line).

**Figure 8.** Seasonal maximum precipitation intensities for 120 minutes of the 7 stations (25% and 75% percentile, median (red)) and the trend 1979–2017 (dashed line).

#### *3.3. River Flow Trends*

The annual mean flow 1981–2016 at the catchment outlet *A* shows a decreasing trend around 1.10 m3s−1, with the maximum mean flow observed in 1998 (3.01 m3s−1) and the minimum mean flow observed in 2016 (0.12 m3s−1) (Figure 9 left).

**Figure 9.** Annual mean (black, left), minimum (blue, **left**) and maximum (red, **right**) flow at Saifenbach and linear trends 1981–2016.

Observed mean annual minimum flows were increasing 1981–2016 around 0.11 m3s−<sup>1</sup> with the smallest recording in 2002 (0.03 m3s−1) and the largest recording in 2014 (0.24 m3s−1) (Figure 9 left). Observed mean annual maximum flows showed no significant trend 1981– 2016 around 31.10 m3s−1, with the largest observation in 1992 (92.14 m3s−1) and the smallest observation in 2015 (5.61 m3s−1) (Figure 9 right).

Mean winter flows 1981–2016 show, as already observed for annual flows, a decreasing trend around 0.51 m3s−1, with the lowest observation in the winter 2016 (0.12 m3s−1) and the largest observation in the winter 1992 (1.99 m3s−1). Minimum winter flows showed no significant trend 1981–2016 around 0.14 m3s−1, with the lowest flow in 2002 (0.03 m3s−1) and the largest minimum in 2014 (0.35 m3s−1). Maximum winter flows also remained constant 1981–2016 at 4.11 m3s−1, with the highest flow recorded in 1992 (34.28 m3s−1) and the lowest maximum in 1984 (0.58 m3s−1) (Figure 10 top).

**Figure 10.** Mean (black), minimum (blue) and maximum (red) flow at Saifenbach and linear trends 1981–2016 for the winter (**top**), spring (2nd from top), summer (3rd from top), and autumn (**bottom**).

Mean spring flows 1981–2016 were decreasing around 0.88 m3s−1, with the lowest mean in 2002 (0.13 m3s−1) and the largest mean recorded in 1994 (5.62 m3s−1). Mean minimum spring flows show no trend at 0.17 m3s−1, with the lowest flow occurring in 2014 (0.33 m3s−1) and the highest minimum observed in 2002 (0.05 m3s−1). Mean maximum spring flows were increasing 1981–2016 around 8.27 m3s−1, with the largest recording in 1994 (42.42 m3s−1) and the smallest recording in 1993 (0.91 m3s−1) (Figure 10 2nd from top).

Mean summer flows 1981–2016 remained constant around 1.64 m3s−1, with the largest summer mean flow observed in 1997 (8.07 m3s−1) and the lowest mean in the summer 2001 (0.19 m3s−1). Summer minima show no trend 1981–2016 at 0.20 m3s−1, with the lowest observation in 2003 (0.04 m3s−1) and the highest in 1986 (1.30 m3s−1). Summer maxima increased 1981–2016 around 26.57 m3s−1, with the largest summer flow in 1992 (92.14 m3s−1) and the lowest maximum in 1984 (0.76 m3s−1) (Figure 10 3rd from top).

Mean autumn flows showed no trend 1981–2016 around 1.08 m3s−1, with the lowest mean recorded in 2001 (0.16 m3s−1) and the largest mean occurring in 1998 (4.35 m3s−1). Autumn minima decreased around 0.19 m3s−1, with the smallest flow recorded in autumn 1992 (0.06 m3s−1) and the largest minimum in 1982 (0.40 m3s−1). Maximum autumn flows remained constant 1981–2016 around 13.04 m3s−1, with the smallest maximum in

autumn 2008 (0.73 m3s−1) and the largest autumn flow in 1998 (60.81 m3s−1) (Figure 10 bottom). A comprehensive summary of the observed runoff trends, including statistical trend properties, is given in Table 4.

#### *3.4. Water Balance and Evapotranspiration*

Particularly in the 1990s, flow measurements at *A* have large gaps preventing the computation of annual flow volumes. Thus, 21 years were available for the assessment of evapotranspiration based on precipitation and catchment runoff (Figure 11). The mean runoff fraction of the water balance 1981–2015 was 55% showing a decreasing trend. It is to be noted though that fewer data were available for the time period 1981–2000 (6 years) than for the period 2001–2015 (15 years). The highest runoff fraction was observed in 2014 with 89%, while the lowest fraction occurred in 2008 with only 30%. In absolute values, the catchment runoff ranged between 131 and 743 mm, with a mean of 338 mm per year.

Based on long-term precipitation and runoff trends, the actual evapotranspiration fraction showed no significant trend 1981–2015, with a mean of 45%, a minimum of 11% in 2014, and a maximum of 70% in 2008. In absolute numbers, actual catchment evapotranspiration 1981–2015 was around 257 mm, with a minimum of 92 mm in 2014 and a maximum of 451 mm in 2008.

**Figure 11.** Annual water balance as fallen precipitation (100%, blue) and runoff fraction (red). The dashed lines mark the long-term trend of the runoff fraction (red) and the actual evapotranspiration fraction (green) 1981–2015. Missing years did not provide sufficient runoff data for a cumulative annual runoff value.

The potential catchment evapotranspiration (PET) was computed using remote sensing data for air temperature, global radiation, and top of the atmosphere solar radiation for the period 1981–2019. Catchment PET rates show an increasing trend 1981–2019 around a mean of 759 mm per year, with a minimum of 711 mm in 1995 and a maximum of 822 mm in 2002 (Figure 12).

**Figure 12.** Annual potential evapotranspiration (PET) based on de Bruin [52] computed with remote sensing data [51] and trend 1981–2019 (dashed line).

#### **4. Discussion**

The mean annual air temperature in the catchment *Pöllau* has been decreasing since 1991, while annual minima remained constant and maxima were increasing. While the development of minima and maxima is a common consequence of ongoing climate change (e.g., [53,54]), the decreasing long-term development of the annual mean temperature in *Pöllau* is less often confirmed by literature (e.g., [20]) as clearly more often rising temperatures are reported (e.g., [15,18,19,55]). For Austria, rising temperatures have also been reported [56], which confirms that climate change impacts at the local or regional scale differ from large-scale assessments. On the other hand, the rising maximum temperatures in the catchment are in line with the APCC report [56]. It is to be noted that the observed time series in *Pöllau* covers approximately 30 years and is thus rather short for temperature change detection. It might, therefore, well be that the analyzed time period coincided with a period where warming in the catchment did not occur (see, e.g., [57]). This assumption is also confirmed by reports and studies addressing climate change in Austria (e.g., [56,58,59]).

The reported climate-change-induced perturbations to precipitation patterns are far more diverse than for air temperature. Increasing [11,23] and decreasing precipitation rates [20,21] were reported as well as areas where no change was detected [19,20,60]. Mean annual precipitation in *Pöllau* remained constant between 1979 and 2019. This observation is confirmed by the Austrian APCC report [56], reporting increased precipitation for the Austrian alpine areas and a decrease for South-East Austria since the beginning of observations. The *Pöllau* catchment falls in between these two areas in the Eastern alpine foothills. The seasonal precipitation analysis indicates a shift towards the summer season, for which an increasing trend was observed. The remaining seasons (spring, autumn, winter) showed no significant trend concerning the fallen precipitation 1979–2019. It is to be noted that the circumstance in which the rain scales are not equipped with heaters allows for the recording of only melted snowfall. Thus, winter and partly spring precipitation trends depend both on snowfall and melting occurring. Seasonal shifts in precipitation have also been reported by earlier studies (e.g., [9,10]), but it is to be noted that especially the climate change induced impact on precipitation shows obvious regional differences [56]. Precipitation intensities for the analyzed durations were increasing for spring and summer. While summer precipitation depth 1979–2019 was also increasing, it remained constant for spring, allowing the assumption of a reduction of events and at the same time, a higher event precipitation. For winter and autumn, no significant trends were detected, as already observed for the precipitation depth in these seasons. Furthermore, for the precipitation intensity trends, it is to be noted that snowfall was only recorded indirectly via melted snow.

Mean river flows at gauge *A* decreased annually as well as for spring and summer, while flow minima increased annually and for autumn, and flow maxima increased for spring and summer. These observations are in line with the APCC report [56]. At the same time, the precipitation depth increased only during the summer season and analyzed precipitation intensities during spring and summer. The rather opposite trends for the precipitation depth and mean river flows indicate that more water is evapotranspirated in the catchment during the warm season, and increasing flow maxima during spring and summer could be due to increasing precipitation intensities for the same seasons.

Based on observed catchment precipitation and runoff, the actual annual catchment evapotranspiration showed no significant trend between 1981 and 2015. It is to be noted though that only river flow at the catchment outlet was used for computation as subsurface flow data were not available. PET rates show an increasing trend 1981–2019 for the catchment. It is to be noted though that these rates were computed based on remote sensing data as local ground climate data were insufficiently complete for PET computation. Furthermore, the selected PET computation approach [52] did not account for wind speed, due to missing data, and may therefore underestimate catchment PET. These observations are

confirmed by several studies reporting similar evapotranspiration trends for the Northern hemisphere [19,29–32] as well as by the APCC report [56].

#### **5. Conclusions**

The presented analyses of hydro-meteorological variables observed in a hydrological research basin in Eastern Austria mostly confirm the results of earlier studies. At the same time, the results confirm the assumption that climate change impacts vary regionally, and large-scale assessments cannot account for site-specific conditions.


**Author Contributions:** Conceptualization, G.K. and D.M.; data curation, G.K.; formal analysis, G.K.; funding acquisition, D.M.; investigation, G.K.; methodology, G.K. and D.C.; project administration, G.K. and D.M.; supervision, D.M.; validation, G.K., D.C. and D.M.; visualization, G.K.; writing original draft, G.K.; writing—review and editing, G.K., D.C. and D.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Restrictions apply to the availability of these data. Data were obtained from the Institute of Urban Water Management, Graz University of Technology and Department 14, Federal State of Styria, and are available from the authors with the permission of the Institute of Urban Water Management, Graz University of Technology and Department 14, Federal State of Styria.

**Acknowledgments:** Roland Fuchs and Christian Rath for operation and maintenance of the test catchment; Department 14 Styria (Robert Schatzl and Josef Quinz) for the support for the test catchment operation. Open Access Funding by the Graz University of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Performance Evaluation and Comparison of Satellite-Derived Rainfall Datasets over the Ziway Lake Basin, Ethiopia**

**Aster Tesfaye Hordofa 1,\*, Olkeba Tolessa Leta 2, Tena Alamirew 3, Nafyad Serre Kawo <sup>4</sup> and Abebe Demissie Chukalla <sup>5</sup>**


**Abstract:** Consistent time series rainfall datasets are important in performing climate trend analyses and agro-hydrological modeling. However, temporally consistent ground-based and long-term observed rainfall data are usually lacking for such analyses, especially in mountainous and developing countries. In the absence of such data, satellite-derived rainfall products, such as the Climate Hazard Infrared Precipitations with Stations (CHIRPS) and Global Precipitation Measurement Integrated Multi-SatellitE Retrieval (GPM-IMERG) can be used. However, as their performance varies from region to region, it is of interest to evaluate the accuracy of satellite-derived rainfall products at the basin scale using ground-based observations. In this study, we evaluated and demonstrated the performance of the three-run GPM-IMERG (early, late, and final) and CHIRPS rainfall datasets against the ground-based observations over the Ziway Lake Basin in Ethiopia. We performed the analysis at monthly and seasonal time scales from 2000 to 2014, using multiple statistical evaluation criteria and graphical methods. While both GPM-IMERG and CHIRPS showed good agreement with ground-observed rainfall data at monthly and seasonal time scales, the CHIRPS products slightly outperformed the GPM-IMERG products. The study thus concluded that CHIRPS or GPM-IMERG rainfall data can be used as a surrogate in the absence of ground-based observed rainfall data for monthly or seasonal agro-hydrological studies.

**Keywords:** CHIRPS; GPM-IMERG; rainfall data scarcity; agro-hydrology; Rift Valley Lake Basin

#### **1. Introduction**

Climate change and variability trend analyses need consistent and long-term time series climate data [1–8] that are required to study the impact of climate change on the agro-hydrological system [9–11]. Such climate studies can benefit from the freely available Global Climate Models (GCMs) outputs such as rainfall data. In addition, complete and long-term rainfall data with high spatial and temporal resolutions are of importance for water resources planning and optimization of crop water productivity especially in waterscarce areas [12–19].

The application of the GCMs rainfall data requires long-term observed-rainfall data for the downscaling and bias correction of coarse resolutions GCMs products into fine resolutions [9,10]. Ground-based rainfall measurement is the most common approach and well recognized as an accurate dataset [20,21]. However, records from the groundbased station are inconsistent over several parts of the world, including Ethiopia [22,23].

**Citation:** Hordofa, A.T.; Leta, O.T.; Alamirew, T.; Kawo, N.S.; Chukalla, A.D. Performance Evaluation and Comparison of Satellite-Derived Rainfall Datasets over the Ziway Lake Basin, Ethiopia. *Climate* **2021**, *9*, 113. https://doi.org/10.3390/cli9070113

Academic Editors: Mohammad Valipour and Sayed M. Bateni

Received: 14 May 2021 Accepted: 28 June 2021 Published: 8 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Furthermore, available weather stations are inadequate and unevenly distributed to capture rainfall spatial heterogeneity, including less accessibility in remote areas [1,24]. This is a prominent problem, especially in developing countries, including the Ziway Lake Basin [25,26].

The advancement and application of remote sensing technologies offer the possibility of using remotely sensed rainfall data in places where ground-based observed rainfall data are not available [24,27–31]. Several satellite-based rainfall products have been developed with promising approaches for obtaining rainfall estimates at regional and global scales, including blending the ground-based observed rainfall data with remotely sensed data [32]. Some of those satellite-based rainfall products include Tropical Precipitation Measuring Mission Multi-Satellite Precipitation Analysis (TMPA) [33], Precipitation Estimation from Remote Sensed Information using Artificial Neural Networks (PERSIAN) [34], Climate Hazards Infrared Precipitation with Stations (CHIRPS) [35], and Global Precipitation Measurement Integrated Multi-SatellitE Retrieval (GPM-IMERG) [36,37].

Globally, several researchers have evaluated the performance of GPM-IMERG rainfall data using ground-based observations or other existing satellite-based rainfall products [28,38–41]. For example, Tong et al. [38] evaluated the monthly performance of the GPM-IMERG rainfall product using gauge observations at both grid and basin scales for the Nanliu River Basin, Beibu Gulf (Southern coast of China). They concluded that the IMERG showed a high accuracy when detecting light rainfall. Anjum et al. [28] demonstrated IMERG-final run rainfall product estimates by comparing it with gauges and TMPA-based real-time data over the northern highlands of Pakistan at annual, monthly, seasonal, and daily time scale. Their study report showed that the IMERG-final run reasonably well performed than the TMPA-based rainfall estimates. Morsy et al. [40] compared TMPA and IMERG rainfall datasets in the arid environment of El-Qaa Plain, Sinai. They concluded that the IMERG data exhibit superior performance than TMPA in all rainfall intensities. Similarly, Kawo et al. [41] evaluated GPM-IMERG early and late run rainfall estimates with ground gauged rainfall at monthly and seasonal time scales over the Lake Hawassa catchment, Ethiopia. They found that both IMERG-early and late run captured the observed rainfall patterns and values during the rainy season than the dry season.

Many studies have also evaluated the performance of CHIRPS and compared it with ground-based observations at different spatial and temporal scales [31,42–50]. For instance, Wu et al. [50] evaluated the performance of the CHIRPS rainfall dataset against groundbased observed rainfall data over the Yunnan Province, China at monthly, annual, and seasonal scales. They found that CHIRPS data performed well in estimating annual and monthly precipitation. Luo et al. [43] evaluated TRMM and CHIRPS rainfall products in the Lower Lancang-Mekong River Basin. They reported that TRMM rainfall products outperformed the CHIRPS rainfall products. Further, Taye et al. [44] evaluated the performance of CHIRPS and Multi-Source Weighted-Ensemble Precipitation (MSWEP) at a monthly time scale over the upper Blue Nile Basin, Ethiopia. They found that CHIRPS better simulated the magnitude of drought than MSWEP in the different elevation zones of the Upper Blue Nile Basin. Goshime et al. [46] conducted a performance evaluation of CHIRPS rainfall product with the gauged rainfall at monthly and daily temporal resolutions over the Lake Ziway Basin, Ethiopia, and concluded that CHIRPS performed better at the monthly time scale. While several studies have been conducted on evaluating the performance of IMERG and CHRIPS, the previous studies have not simultaneously evaluated and compared the performance of the three IMERG runs (early, late, and final) and CHIRPS at different time scales (monthly and seasonal). Therefore, evaluating and comparing the performance of the recently available different rainfall products at two-time scales is of interest for in-depth and better understanding of their performance and appropriately choosing them as a surrogate when ground-based rainfall observations are lacking. Such studies might also help to identify at what time resolution the satellite-based rainfall estimates can appropriately be used as they play a key role in simulating long-term agro-hydrological modeling and in forecasting changes in freshwater supply and agricultural crop yields [51,52]. Thus, the

objectives of this study were to evaluate the accuracy of the satellite-based areal rainfall data over the Ziway Lake Basin at different time scales. We evaluated and compared the CHIRPS and GPM-IMERG of early, late, and final runs with the ground-based observed rainfall data from 12 gauging stations. The evaluation was performed at monthly and seasonal time scales from 2000 to 2014. This study might be useful for the alternative application of remotely sensed precipitation products in simulating the agro-hydrological modeling and climate change trend assessment of the Ziway Lake Basin and elsewhere with similar agro-hydrological conditions, in the Central Rift Valley Lake Basin of Ethiopia.

#### **2. Data and Methods**

#### *2.1. Study Area Description*

Lake Ziway Basin (LZB) is located between 38◦00 −39◦30 East longitude and 7◦00 −8◦30 North latitude in the Adami Tullu-Jiddo Kombolcha Woreda of the East Shewa Zone, Oromia region, Ethiopia. The basin is about 150 km south of the capital city, Addis Ababa. The town of Ziway (recently named Batu) is situated on the lake's western shore. The altitude of Lake Ziway is approximately 1636 m above mean sea level (amsl), with a maximum water depth of 4 m, a total basin area of about 7300 km<sup>2</sup> (Figure 1) and a lake volume of 1.5 million cubic meters [53]. The majority of the basin is characterized by low to moderately undulating topography but bounded by a steep slope and abrupt faults in the eastern and southeastern escarpments, ranging from 4200 to 1600 m (Figure 1). Lake Ziway Basin experiences the monsoon agro-climate zone characteristics. The rainfall patterns are generally affected by the annual oscillation of the inter-tropical convergence zone that forms wet summer from June to September [54]. The mean annual rainfall of the basin spatially varies from 500 to 1150 mm, with a noticeable temporal variation at a monthly time scale. The mean annual temperature ranges from approximately 15 ◦C for the highlands to 25 ◦C close to the lake.

**Figure 1.** A map of the Ziway Lake Basin, including elevation, rivers, rainfall stations, and Lake Ziway itself.

#### *2.2. Data*

#### 2.2.1. Ground Observed Data

In this study, the monthly and seasonal rainfall ground-based observed data from 2000 to 2014 were used as a point of reference for evaluating the CHIRPS and GPM-IMERG. We obtained the data from the Ethiopian National Meteorological Agency (NMA). We

originally obtained nineteen climate stations distributed over the Ziway Lake basin with different elevation. However, after performing quality and checking consistency of the data, we selected 12 stations that had good quality and consistent temporal coverage (Table 1). Then, we applied Thiessen polygon method in order to calculate the areal weighted rainfall values of the Ziway Lake Basin (ZLB) from the 12 selected stations. Such approach accounts for the areal coverage of each rain gauge station, the spatial distribution and variability of rainfall for the basin [55]. The areal coverage (Thiessen polygon) of the 12 stations is shown in Figure 2.


**Table 1.** List of the twelve rainfall stations over the Ziway Lake Basin.

**Figure 2.** Thiessen Polygon network of the Ziway Lake Basin.

#### 2.2.2. Satellite Precipitation Products

In this study, we considered and evaluated two Satellite Precipitation Products (SPPs). These are CHIRPS and GPM-IMERG.

#### CHIRPS Database

CHIRPS was launched in early 2014 by the Climate Hazards Group at the University of California, Santa Barbara (UCSB). The CHIRPS precipitation dataset globally covers 50◦ S−50◦ N with a horizontal resolution of 0.05◦ for both daily and monthly time scales. CHIRPS datasets were originally developed to support the United States Agency for

International Development Famine Early Warning Systems Network (FEWS NET) [35] and African Rainfall Climatology [55,56]. Nowadays, the CHIRPS dataset is available in two sets of spatial resolutions i.e., 0.25◦ × 0.25◦ and 0.05◦ × 0.05◦ from 1981 to the present.

The CHIRPS dataset is developed based on a blend of three data sources [35]: (i) the Climate Hazards Precipitation Climatology (CHPclim) [57], a global precipitation climatology at 0.05◦ latitude and longitude resolution (estimated for each month based on station data, averaged satellite observations, elevation, latitude and longitude) [35,58]; (ii) quasi-global geostationary Thermal Infrared Radiation (TIR) satellite observations, TMPA 3B42 product [33], and (iii) atmospheric model precipitation fields from the National Oceanic and Atmospheric Administration (NOAA) Climate Forecast System (CFS) version 2.0 [59].

According Funk et al. [35], the CHIRPS algorithm encompasses four development processes: (i) a pentad (5 day) rainfall estimate, which is generated from the three-hourly quasi-global geostationary TIR data of Climate Prediction Center (CPC) and the National Climatic Data Center; (ii) a TMPA-3B42 rainfall product, which is used to calibrate the IR pentad estimate; (iii) the calibrated IR pentad product is then multiplied with the Climate Hazards Precipitation Climatology and subsequently divided by the longterm mean to produce the Climate Hazards Group (CHG) IR Precipitation (CHIRP) data; (iv) the pentadal CHIRP values are disaggregated to daily precipitation estimates based on the daily NOAA Climate Forecast System (CFS) fields rescaled to 0.05◦ resolution. Finally, CHIRPS is produced through blending the rainfall stations with the CHIRP data sets and using a modified inverse distance-weighted algorithm [35].

The CHIRPS datasets include rainfall information from a large number of gauges, which is about 1200 stations globally. It should be mentioned that a relatively large number of rain gauge stations were used in East Africa [35]. More than 50 rain gauge stations from the Ethiopian NMA were blended with the CHIRPS products for up-to-date evaluations of the rainfall conditions throughout the major growing seasons of the country. The 50 stations are updated every 10 days [60] and used to correct the CHIRPS datasets [35,49,61] Detailed information regarding the CHIRPS rainfall products was provided in Funk et al. [35]. In this study, we used a higher resolution CHIRPS dataset with a spatial resolution of 0.05◦ × 0.05◦ and a daily time scale, which was freely downloaded from (ftp://ftp.chg. ucsb.edu/pub/org/chg/products/CHIRPS-2.0/).

#### IMERG Database

The GPM-IMERG algorithm combines information from the GPM satellite group to estimate precipitation over the majority of the Earth's surface. The GPM-IMERG was launched by the National Aeronautics and Space Administration (NASA) and the Japan Aeronautics and Exploration Agency (JAXA) in 2014 [62]. This algorithm is particularly valuable over the majority of the Earth's surface that lacks precipitation-measuring instruments on the ground. In the latest release of IMERG (Version 06; V06), the algorithm fuses the early precipitation estimates based on the TRMM satellite (2000−2014) with more recent precipitation estimates collected during the operation based on the GPM satellite (2014–2021). The three gridded products are commonly used for scientific research and operational purposes. There are three different daily IMERG products, which include IMERG Day 1 Early Run (near real-time with a latency of 6 h), IMERG Day 1 Late Run (reprocessed near real-time with a latency of 18 h), and IMERG Day 1 Final Run (gauged-adjusted with a latency of four months) products. In this study, we used the three IMERG products (IMERG-early IMERG-late and IMERG-final run products, with a fine spatial resolution (0.1◦ × 0.1◦), a high temporal resolution (30 min), and a spatial coverage from 60◦ S to 60◦ N, which was freely downloaded from (https://giovanni.gsfc.nasa.gov/giovanni/ (accessed on 4 February 2021)).

#### *2.3. Performance Evaluation Criteria*

To identify the best datasets in the study area, we evaluated the performance of CHIRPS and three IMERG (early, late, and final) products against the ground-based rainfall data. We evaluated the monthly and seasonal time scale. We obtained monthly and seasonal rainfall by adding up the daily values on a monthly and seasonal basis in Microsoft Excel 2019 [63], Jupyter Notebook and ArcMap used to visualize data. In Ethiopia, the climate varies mostly with altitude. The lowland areas have hot and arid climatic conditions while plateau areas experience a cold climate, and the season category does not constant over the regions [64,65]. Therefore, in this study, we characterized the performance of CHIRPS and IMERG rainfall datasets for the four seasons of the ZLB. These include *Kiremt* (summer; from June to August), *Tseday* (spring; from September to November), *Bega* (winter; from December to February), and *Belg* (Autumn; from March to May). Then, we evaluated the temporal variations of rainfall for each product.

We consistently used four statistical metrics that include Percent Bias (PBIAS), Root Mean Square Error (RMSE), Nash–Sutcliffe Efficiency (NSE), and Pearson linear Correlation Coefficient (r) to quantitatively compare the performance of the CHIRPS and the three GPM-IMERG rainfall products. PBIAS describes the systematic bias of the CHIRPS and IMERG products. Positive values of PBIAS indicate an overestimation of the rainfall quantity, whereas negative values show an underestimation of the rainfall quantity [28,66,67]. RMSE measures the absolute error magnitude of the CHIRPS and IMERG products, with the smaller the RMSE value, the closer the CHIRPS and IMERG measurements to the groundobserved rainfall. NSE is a normalized statistic that determines the relative magnitude of the residual variance compared to the measured data variance. NSE values range between −∞ and 1, with value 1 indicating a perfect fit between the satellite-based and observed rainfall [42,68]. The degree of linear correlation between the CHIRPS and IMERG and the ground-based rainfall evaluated with r values ranging from −1 to 1 r value of 0 indicates no correlation between the CHIRPS and IMERG products and the observed rainfall. On the other hand, r values of 1 and −1 show perfect positive and negative correlations, respectively [69,70], as summarized in (Table 2). In addition to statistical metrics, we used graph for comparison of SPPs and observed rainfall.


**Table 2.** List of the statistical metrics, used for the evaluation of satellite rainfall products.

where: PSi is rainfall from satellite and PGi the observed rainfall at ith time step (daily, weekly, monthly, or seasonal) with N pairs of data, PGmean and PSmean are mean observed rainfall and mean satellite rainfall, respectively.

#### **3. Results and Discussion**

#### *3.1. Spatial Rainfall Pattern Evaluation*

The Ziway Lake Basin seasonal average rainfall distribution of the CHIRPS and IMERG map was compared visually from the 2000–2014 period. Figure 3 shows the seasonal average rainfall distribution for the main rainy (summer) and dry (winter) seasons.

In summer (Figure 3a,b), both CHIRPS and IMERG show that the western part of the basin, which is the eastern highlands of Gurage Zone, receives more rainfall than the eastern part of the basin, which is the western highlands of the Arsi Zone. The spatial rainfall distribution of both CHRIPS and IMERG is consistent with ground-observed rainfall [64]. During the winter season (DJF), a similar rainfall pattern was observed in the western and eastern parts of the basin (Figure 3c,d). Up to 105 mm of rainfall amount is received for the eastern and western part of the basin whereas the central and southern part of the basin receives rainfall up to 45 mm. Overall, both CHRIPS and IMERG showed a decreasing rainfall pattern towards the center i.e., from west to the central part of Ziway Lake Basin (lowland). According to Hailesilassie et al. [64], the observed rainfall is mainly concentrated in the southern and western parts of the basin, while the eastern and central rift valley (low land areas) where the lake is located generally experience low rainfall amounts. CHIRPS relatively well captured that pattern when compared to IMERG, which is probably due to its high spatial resolution and blending of more stations' data [47].

**Figure 3.** Spatial distribution of main rainy and dry season rainfall (**a**,**c**) for CHIRPS, (**b**,**d**) for IMERG for the period 2000–2014.

#### *3.2. Monthly Rainfall Evaluation*

Comparison of the CHIRPS and IMERG (early, late, and final run) monthly rainfall data showed a good performance over the Ziway Lake Basin. CHIRPS rainfall generally showed a stronger correlation with the observed rainfall when compared to the three-run IMERG's rainfall (Table 3). The Correlation Coefficient between the early, late, and final IMERG run rainfall and the observed rainfall was high i.e., 0.93, 0.92, and 0.85, respectively. Compared with all IMERG (early, late, and final) products, CHIRPS products showed the highest Correlation Coefficient (0.96) and low Percent Bias (2.22%). In comparison with the IMERG products, the monthly CHIRPS product relatively better represented the groundobserved rainfall values over ZLB with relatively higher r and NSE; and lower RMSE and RBIAS. This is consistent with the previous studies of that confirmed the applicability of CHIRPS precipitation datasets at a monthly time scale in ground-observed data-scarce regions [22,31,46,49,64].


**Table 3.** Monthly statistical performance evaluation satellite rainfall products for the Ziway Lake Basin.

Figure 4a shows the monthly rainfall values while Figure 4b, c shows the cumulative and scatter values, respectively. The CHIRPS and IMERG-L rainfall product showed the best performance to capture the temporal pattern of monthly rainfall. However, both IMERG-E and IMERG-F products did not well capture the temporal variability of observed rainfall over the study area, indicating that both somehow overestimated the observed rainfall values. As visualized from the cumulative rainfall (Figure 4b), the CHIRPS and IMERG-L captured the monthly cumulative observed rainfall values. The IMERG-E and IMERG-F run smoothly captures the temporal cumulative observed rainfall compared to the CHIRPS and IMERG-L product. As the scatter plot (Figure 4c) indicated, the monthly CHIRPS and IMERG-L rainfall values are close to the monthly observed rainfall values. The CHIRPS data showed capability to represent the monthly maximum observed values compared to all the IMERG's runs. IMERG-L data generally outperformed the IMERG-E and IMERG-F data (Table 3).

**Figure 4.** *Cont.*

**Figure 4.** Monthly areal rainfall (**a**), cumulative rainfall depths (**b**), the correlation between monthly satellite-derived (CHIRPS and IMERG-(early, late, and final) run) and observed rainfall (**c**) from (2000–2014) over the Ziway Lake Basin.

#### *3.3. Seasonal Rainfall Evaluation*

Figure 5 shows statistical metrics used for seasonal rainfall evaluation of the SPPs versus the ground stations. There were some slight differences between these products on r, RMSE, NSE, and PBIAS (Figure 5a–d). The figure shows that the CHIRPS, the IMERG-E, IMERG-L, and IMERG-F performed well. Moreover, the IMERG-E, IMERG-L, and IMERG-F performance indicated a better relationship during the summer season with an r and NSE values of (0.96 and 0.9 and (0.95 and 0.96), respectively, whereas CHIRPS well-performed with a high r value of 0.92 and low bias error (−2.6) (Figure 5a–d). The three IMERG runs underestimated the summer rainfall by −2.9% to −10%, while CHIRPS underestimated the summer season rainfall by −12% (Figure 5d). All IMERG runs overestimated observed rainfall by 4% to 9.7% in the winter season, whereas CHIRPS underestimated the observed values by −2.6% (Figure 5d). When compared to IMERG runs, CHIRPS achieved higher correlations with observed rainfall during spring, winter, and autumn seasons with r values of 0.93, 0.97, and 0.93 (Figure 5a), respectively. The RMSE values indicated that the CHIRPS data relatively had a small value compared to all IMERG runs, especially during the winter and autumn seasons (Figure 5b). During the spring season, the three IMERG runs had the same r values (0.92) and CHIRPS had (0.93) (Figure 5a).

**Figure 5.** Seasonal performance evaluation indices of CHIRPS, IMERG-E, IMERG-L, and IMERG-F run: Correlation Coefficient (**a**), Root Mean Square Error (**b**), Nash–Sutcliffe Efficiency (**c**), and Percent Bias (**d**) for the period 2000–2014.

A number of previous studies reported the good performance of SPPs at monthly time scales [25,28,41,46,50,68–70]. In general, CHIRPS showed slightly better performance than the other three IMERG runs for monthly and seasonal time scales. Previous studies have already confirmed the superiority of CHIRPS than IMERG runs for different parts of the world [40,71–73], including Ethiopia [31,47,49]. For example, Wedajo et al. [47] reported better rainfall estimation by CHIRPS compared with IMERG and TAMSAT3 and 3B42/3 products for the Dhidhessa River Basin, Ethiopia. Dinku et al. [49] reported better rainfall estimation capability of CHIRPS for east Africa compared to the African Rainfall Climatology version 2 (ARC2) and TAMSAT3 products. The better performance of CHIRPS has been attributed to the capability of the algorithm to integrate satellite, rain gauges, and reanalysis products, combined with its higher spatial and temporal resolutions than IMERG products [35].

Overall, the statistical evaluation results indicate that both CHIRPS and IMERG are capable of estimating and detecting observed monthly and seasonal rainfall values of the ZLB. Therefore, the monthly and seasonal CHIRPS and IMERG-F data are a reliable source for simulating monthly and seasonal agro-hydrological processes, estimating the seasonal crop water requirement, and accounting the stocks and fluxes of water in the Ziway Lake Basin.

#### **4. Conclusions**

In this study, we evaluated and compared the performance of IMERG and CHIRPS rainfall products against ground-observed rainfall data over the Ziway Lake Basin. The analyses covered the period from 2000 to 2014 at monthly and seasonal time scales. We used four statistical evaluation parameters: Correlation Coefficient, Nash–Sutcliffe Efficiency, Percent Bias, and Root Mean Square Error. The two rainfall products performed well for both monthly and seasonal time scales. Overall, while the CHIRPS's rainfall datasets showed slightly better performance over the IMERG's datasets, both datasets can be used at a monthly or coarser temporal resolution when ground-based rainfall data are not available. This can greatly contribute to continuous spatiotemporal monitoring of drought and helping the water managers and agricultural planners implementing mitigation measures and improving the livelihood of the stakeholders in the basin.

The follow up research should focus on the evaluation and comparison of the grid point satellite dataset with interposed ground station data, considering point to point performance evaluation at daily time basis. Future evaluation studies should also include the Climate Hazards Group Infrared Precipitation (CHIRP) satellite-only product.

**Author Contributions:** Conceptualization, A.T.H.; data curation and analysis, A.T.H. and N.S.K.; Formula analysis and Methodology, A.T.H.; writing–original draft, A.T.H.; writing–review, A.T.H.; editing, A.T.H., O.T.L. and A.D.C.; supervision O.T.L., T.A. and A.D.C. All authors have read and agreed to published version of the manuscript.

**Funding:** This work was financially supported by the Africa Center of Excellence for Water Management, Addis Ababa University.

**Acknowledgments:** The authors would like to thanks the National Meteorological Agency (NMA) of Ethiopia for providing the weather data. We would also like to express our sincere gratitude to the Africa Centre of Excellence for Water Management, Addis Ababa University for the support to conduct this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


# *Article* **Modeling the Soil Response to Rainstorms after Wildfire and Prescribed Fire in Mediterranean Forests**

**Manuel Esteban Lucas-Borja 1, Giuseppe Bombino 2, Bruno Gianmarco Carrà 2, Daniela D'Agostino 2, Pietro Denisi 2, Antonino Labate 2, Pedro Antonio Plaza-Alvarez <sup>1</sup> and Demetrio Antonio Zema 2,\***


Received: 17 November 2020; Accepted: 15 December 2020; Published: 17 December 2020

**Abstract:** The use of the Soil Conservation Service-curve number (SCS-CN) model for runoff predictions after rainstorms in fire-affected forests in the Mediterranean climate is quite scarce and limited to the watershed scale. To validate the applicability of this model in this environment, this study has evaluated the runoff prediction capacity of the SCS-CN model after storms at the plot scale in two pine forests of Central-Eastern Spain, affected by wildfire (with or without straw mulching) or prescribed fire and in unburned soils. The model performance has been compared to the predictions of linear regression equations between rainfall depth and runoff volume. The runoff volume was simulated with reliability by the linear regression only for the unburned soil (coefficient of Nash and Sutcliffe E = 0.73–0.89). Conversely, the SCS-CN model was more accurate for burned soils (E = 0.81–0.97), also when mulching was applied (E = 0.96). The performance of this model was very satisfactory in predicting the maximum runoff. Very low values of CNs and initial abstraction were required to predict the particular hydrology of the experimental areas. Moreover, the post-fire hydrological "window-of-disturbance" could be reproduced only by increasing the CN for the storms immediately after the wildfire. This study indicates that, in Mediterranean forests subject to the fire risk, the simple linear equations are feasible to predict runoff after low-intensity storms, while the SCS-CN model is advisable when runoff predictions are needed to control the flooding risk.

**Keywords:** hydrological models; rainfall; surface runoff; linear regression models; curve number; SCS.CN model; mulching; wildfire; prescribed fire

#### **1. Introduction**

The importance of forests in the climatic context at the planetary scale is well known, since forests produce oxygen and store carbon, regulate water and energy fluxes, support biodiversity and provide other fundamental ecosystem services [1–3]. However, the fundamental role of forests is threatened by some natural and anthropogenic agents, such as the extreme weather events and fire, with a long history of influence on forest ecosystems [4]. Extreme weather events (e.g., heavy storms and severe drought) are more and more intensified by climate change trends and occur in all regions of the world [5,6]. The fire effects extend to several components of forests, such as soil, vegetation, air and surface water [7]. With regard to the effects on surface water and soil, fire strongly alters the hydrological response of recently burnt areas, increasing by many folds the soil's aptitude to generate runoff and erosion compared the unburned forest areas [8,9]. High runoff and erosion rates in

forests increase the flooding risk, debris flow occurrence and water quality alteration in downstream areas, with possible loss of human lives and heavy damage to infrastructures and environment [10,11]. For instance, with regard to the Mediterranean climate, exceptionally high erosion rates (up to 100 tons per ha) have been reported immediately after wildfire by Menendez-Duarte et al. [12] in the Iberian peninsula, while Lopez-Batalla et al. [13] measured increases in flood runoff by 30% and in peak discharges by 120% in the same environments. These studies together with a large body of literature (see several examples in the milestone review of Shakesby [14]) clearly demonstrate the need to control and mitigate the hydrological response of forest soils after the wildfire adopting prediction tools and post-fire management actions.

The hydrological processes in forests are influenced by several factors, among which fire severity is important [15]. In other words, the more severe the fire, the more susceptible the soil to increases in runoff and erosion and worsening in water quality changes in the downstream ecosystems [7]. For instance, Lucas-Borja et al. [7] found increases in runoff and erosion in soil burned by wildfire by about 20% and even 200%, respectively, compared to unburned soils in Central-Eastern Spain.

The hydrology of burned forests is very complex, since it is the product of several factors, such as climate and edapho-climatic conditions, fire severity, soil, vegetation, morphology and land management after fire [15–19]. The needs to predict the impacts of fire on runoff and erosion and to implement the most effective actions for rehabilitation of fire-affected areas have increased the demand for hydrological models [16,20]. This demand is particularly important for forest managers working in the Mediterranean Basin, which is characterized by dry and hot summers followed by frequent and high-intensity rains immediately after the wildfire season [15,21]. In Mediterranean areas, increases in wildfire frequency and burned areas are expected under the forecasted climate scenarios [22,23].

In Mediterranean forests, the literature reports applications and verifications of hydrological models with different nature and complexity: from the simplest empirical models (such as the Soil Conservation Service (SCS)-curve number model to predict runoff, the universal soil loss equation, USLE, to simulate soil erosion) through the semiempirical models (e.g., the Morgan–Morgan–Finney model, MMF) until the most complex physically-based models (for instance, the Water Erosion Prediction Project, WEPP) or even the artificial neural networks [24,25]. Nonetheless, the empirical models are sometime more commonly used compared to the more complex models, mainly in data-poor environments (that is, in those situations with limited availability of parameter inputs) and for quick identification of sources of water, sediments and pollutants in forests [24–26].

Accurate predictions of surface runoff are fundamental to achieve reliable estimations of erosion rates and water quality parameters using hydrological models [27], particularly in forests subject to climate change and fire, since the latter factors play a large influence of the soil's hydrological response [28,29]. The Soil Conservation Service (SCS)-curve number (CN) model (hereinafter "SCS-CN model") is one of the most common methods for estimating the runoff volume generated by a rainstorm [30,31]. The popularity of this method is due to its simplicity, ease of use, widespread acceptance and large availability of input data [32]. Moreover, the SCS-CN model takes into account most of the factors that influence runoff generation, such as soil type, land use, hydrologic condition and antecedent moisture of the soil [33]. The model has also been incorporated as a hydrological submodel in several distributed rainfall–runoff models at the watershed scale (e.g., AnnAGNPS—annual agricultural non-point source pollution model, SWAT—soil and water assessment tool model and HEC-HMS—hydrologic engineering center-hydrologic modeling system model), supporting its robustness and popularity [33,34]). To date, there is no any alternative model that offers as many advantages as the SCS-CN model, which therefore is still commonly used in the large majority of environments and climatic conditions [35].

However, various studies conducted throughout the world on the applicability of the SCS-CN method have suggested a need for further improvement or overhauling of the model [32,36], since in some environments the method provides unsatisfactory predictions, particularly when the soil's hydrological response does not follow the runoff generation mechanism by saturation-excess. Moreover, in spite of the large number of requisites, the model has been surprisingly little used for hydrological predictions in fire-affected areas, and the CN values are not completely known in burned areas [37]. The hydrological research has been mainly carried out at the watershed scale, where post-fire runoff has been predicted using the SCS-CN method incorporated in watershed-scale models (e.g., WILDCAT4 Flow Model [38] and FIRE HYDRO [39]). For instance, the SCS-CN model was used by Candela et al. [40], who analyzed the flood frequency curves for pre- and post-fire conditions, showing an increase in the average curve numbers and a decrease in the catchment time lag. Increases by 25 units in post-fire CNs were estimated by Soulis [37] in a small Greek watershed using pre-fire and post-fire rainfall–runoff datasets. A daily-constant CN in the SWAT model was used by Nunes et al. [41] to simulate the effects of soil water repellency on runoff from burnt hillslopes in a Mediterranean forest throughout three years after fire. It is therefore evident that the modeling experiences are scarce at the plot scale. At this scale, modeling of soil hydrology is less complex compared to the watershed scale, where the hydrological response to Mediterranean storms is further complicated by a combination and overlaying of several hydrological processes (e.g., water routing in the channel network, ponding and uneven soil properties) other than surface runoff generation. Furthermore, the studies about the hydrological effects of post-fire management on runoff in forests using the SCS-CN model are scarce.

Therefore, there is a need of further studies that must evaluate the runoff prediction capacity of the SCS-CN model in forest hillslopes or plots affected by fire of different severity—a fire parameter referring to the effects of wildfire on plant communities—in comparison with simpler models that estimate runoff directly from precipitation, such as the linear regression equations. In other words, is the SCS-CN model accurate to predict surface runoff in Mediterranean burned forests? Is it able to simulate post-fire hydrology with or without rehabilitation management actions (such as straw mulching) after a wildfire? When is its use convenient compared to a simpler linear regression between rainfall and runoff?

This study aims to reply to these questions, evaluating the hydrological performance of the SCS-CN model in two pine forests of Central-Eastern Spain affected by a wildfire and a prescribed fire, respectively, having different fire severity. More specifically, observations of rainfall–runoff patterns collected throughout one year in undisturbed soils (assumed as control) and in plots subject to prescribed fire/wildfire (the latter with or without a mulching treatment) are compared with the corresponding predictions of the SCS-CN model and linear rainfall–runoff regressions. The outcomes of this study help land managers to adopt strategies to control the hydrological effects of fire in Mediterranean forests.

#### **2. Materials and Methods**

#### *2.1. Study Areas*

Two experimental areas were selected in pine forests of the Province of Albacete, Castilla—La Mancha Region, Central Eastern Spain. The first area (Sierra de las Quebradas, municipality of Liétor) was affected by a wildfire in July 2016. The second forest (municipality of Lezuza) was subjected to a prescribed fire in March 2016, to reduce fuel loading and thus the potential risk and severity of subsequent fires (Figure 1). Prior to wildfire, the soil cover of the forest was mainly composed of plants, litter and stones with variable composition.

Both study areas have a semiarid Mediterranean climate, BSk according to the Köppen–Geiger classification [42]. The average annual rainfall and medium annual temperature are 282 (Liétor) and 450 (Lezuza) mm and 13.5 (Liétor) and 16 (Lezuza) ◦C, respectively (Spanish National Meteorological Agency, 1950–2016). According to historical data (1990–2014) provided by the Spanish Meteorological Agency, the maximum precipitation is concentrated in October (44.5 mm) and the minimum in May (39.6 mm); from June to September a hot and dry period (air relative humidity below 50%) occurs. The mean minimum temperature of the coldest month is −0.9 ◦C and the mean maximum temperature of the hottest month is nearly 32 ◦C.

**Figure 1.** Location and layout of forest plots subject to prescribed fire and wildfire and monitored for hydrological observations (Lezuza and Liétor, Castilla La Mancha, Spain). Geographic coordinates and map source: Lezuza X: 557588E, Y: 4306475N; Liétor: X: 600081 E, Y: 4262798 N (unburned area); X: 598358 E, Y: 4264032 N (burned area); Google Earth, last access on 6/15/2019).

#### 2.1.1. Wildfire-Affected Forest (Liétor)

The experimental area of Liétor is located at an elevation between 520 and 770 m a.s.l. with W-SW and N aspect and mean slope of 15–20% (Table 1). Soils are classified as Inceptisols and Aridisols with sandy-loam texture [43].

The wildfire burned about 830 ha of the forestland (mainly *Pinus halepensis* Mill). The mean tree density of this forest was between 500 and 650 trees per hectare with height between 7 and 14 m. *Rosmarinus o*ffi*cinalis* L., *Brachypodium retusum* (Pers.) Beauv., *Cistus clusii* Dunal, *Lavandula latifolia* Medik., *Thymus vulgaris* L., *Helichrysum stoechas* (L.) Moench, *Macrochloa tenacissima* (L.) Kunth, *Quercus coccifera* L. and *Plantago albicans* L. are the shrub or herbaceous species of the forest. The wildfire, classified as high-severity fire by the local forest managers according to the methodology proposed by Vega et al. [44], determined a tree mortality of 100% (Table 1). Forest floor was about 3–5 cm deep. This forest floor, as happened also in Lezuza, was blanketed with decaying *Pinus halepensis* M. needles and twigs and other wood debris such as cones or branches coming from trees. In both sites, the forest floor was blanketed with decaying *Pinus halepensis* M. needles and twigs and other wood debris such as cones or branches coming from trees in both sites, Lezuza and Liétor. More information related to this suggestion is provided in Table 1.

In September 2016 (two months after the wildfire), a mulching treatment was carried out in some areas of the burned site. Barley straw was spread manually on soil at a depth of 3 cm and a rate of 0.2 kg m−<sup>2</sup> (dry weight), following the dose proposed by different authors for forests of Northern Spain, to achieve a burned soil cover over 80% [45].


**1.**Maincharacteristicsofforestsandplotssubjecttoprescribedfireandwildfire(LezuzaandLiétor,CastillaLaMancha,

#### 2.1.2. Forest Subjected to Prescribed Fire (Lezuza)

The forest area was selected in a relatively hilly area at an elevation from about 1010 to 1040 m a.s.l., with a 50-year old mixed plantation of *Pinus halepensis* and *Pinus pinaster*. The mean slope is around of 15% and the aspect is N-NE. Soils are classified as Alfisols with Xeralf Rhodoxeralf horizon with clay texture [43] (Table 1).

The tree density of this forest area was about 500 trees per hectare with a mean height of 6.40 m. The understory was dominated by *Quercus faginea* Lam. L., *Quercus ilex* subsp. ballota, *Quercus coccifera* L., *Juniperus oxycedrus*, *Brachypodium retusum* P. and *Thymus* sp. (Table 1). Forest floor depth was about 5–7 cm.

#### *2.2. Description of Experimental Plots and Measurement of Runo*ff *Volume*

In the two selected forests, experimental plots were installed in unburned (control) and burned areas. Plots were randomly distributed in the experimental site in areas with the same morphological and ecological characteristics to ensure comparability.

More specifically, in the wildfire-affected forest of Liétor, eighteen rectangular plots (20-m long × 10-m wide) were installed at a distance between 200 and 500 m. Six plots were located in the forest outside of the burned site and assumed as a control. Twelve plots were instead located in the burned area, of which six were not treated, while six plots were mulched. In the forest subjected to the prescribed fire in Lezuza, twelve plots (4-m long and 2-m wide) were isolated, of which six were located in the unburned area and the other six in the burned site. The prescribed fire was carried out under controlled air conditions in the forests (wind speed of 14 km/h, air temperature of 14 ◦C and relative humidity of 63%), which are reference values for applying the prescribed fire as a forest protection measure. The upper and side borders of all plots in both areas were hydraulically isolated by geotextile fabric pounded into the ground, to prevent external water inputs. At the plot bottom, runoff was collected using a metal fence conveying the water into a pipe, which discharged to a plastic tank of 25 (Liétor) or 50 (Lezuza) liters. In these plots, immediately the runoff volumes were measured after each rainfall event throughout an observation period of about one year (Table 2).


**Table 2.** Main characteristics of precipitation events in plots subject to prescribed fire and wildfire (Lezuza and Liétor, Castilla La Mancha, Spain).

A weather station (WatchDog 2000 Series model) with a tipping bucket rain gauge was placed in each study area to measure total daily precipitation, storm duration, air temperature and rain intensity during the study period. In the hourly rainfall series of the experimental database, two consecutive events were considered separate, if no rainfall was recorded for 6 h or more [46,47]. The mean rainfall intensity was the total rainfall divided by the storm duration. The main variables characterizing the observed events are reported in Table 2.

#### *2.3. Outlines on the SCS-CN Model*

This model, proposed by the Soil Conservation Service of the United States Department of Agriculture in 1972, hypothesizes that:

$$\frac{V}{P\_n} = \frac{W}{S} \tag{1}$$

being:


*Pn* is the difference between the total precipitation (*P*) and the initial losses (*Ia*, as the water storage in the soil dips, interception, infiltration and evapotranspiration) prior to surface runoff. *Ia* is assumed to be proportional to S through a coefficient λ:

$$I\_d = \lambda S \tag{2}$$

*S* is given by:

$$S = 25.4 \cdot \left(\frac{1000}{CN} - 10\right) \tag{3}$$

when the parameter CN is the so called "curve number". The CN can be considered as the soil's aptitude to produce runoff and is a function of the hydrological properties and conditions of soil, and land use. The CN varies between 0 and 100 (0 means that the soil does not produces runoff, 100 means that all the precipitation turns into surface runoff and then the hydrological losses are zero).

According to this model, the runoff volume *V* is:

$$V = \frac{\left(P - \lambda S\right)^2}{P + (1 - \lambda)S} \tag{4}$$

To estimate CN in agroforest areas, the soil hydrological class, vegetation cover, hydrological condition (good, medium and poor) and cultivation practice and the antecedent moisture condition (AMC) of the soil must be determined.

The soil hydrological class (A to D) is related to the soil's capability to produce runoff, on its turn due to the soil infiltration capacity. A low runoff production capability corresponds to the A soil hydrological class, while the highest runoff capability is typical of less permeable soils D.

The actual AMC of the soil subject to a rainfall/runoff event was estimated as a function of the total height of precipitation in the five days before the event in the two different conditions of crop dormancy or the growing season. On this regard, three AMCs are identified:


The SCS guidelines make the CN values available in tables for a soil of given hydrological class and condition, vegetal cover, cultivation practice and average AMC (AMCII). The values of CNs related to AMCI (*CNI*) or AMCIII (*CNIII*) can be calculated through the following equations:

$$\text{CN}\_{I} = \frac{4.2 \text{CN}\_{II}}{10 - 0.058 \text{CN}\_{II}} \tag{5}$$

$$\text{CN}\_{III} = \frac{23 \text{CN}\_{II}}{10 + 0.13 \text{CN}\_{II}} \tag{6}$$

The parameter AMC takes into account the influence of the soil water content on the hydrological response of the soil to the rainstorm and distinguishes "dry" (AMCI), "average" (AMCII) and "wet" (AMCIII) conditions depending on the total rainfall height of the five days before the event.

#### *2.4. Model Implementation*

#### 2.4.1. Linear Regression between Rainfall and Runoff

A linear regression model was established between the surface runoff volume (dependent variable) and the rainfall height (independent variable) for each event, as follows:

$$V = aP\tag{7}$$

where:


The intercept of this linear equation was forced to zero, in order to avoid runoff without any precipitation.

#### 2.4.2. SCS-CN Model

The SCS-CN model was first applied considering the "default" input parameters, that is, the values of λ and CN derived from the SCS guidelines for woods (control plots) or pasture (burned plots) for the soil hydrologic group A of the experimental soils and AMC "I" for all the modeled events (since no or vey low precipitation was recorded in the antecedent five days). However, the runoff prediction capacity was totally unsatisfactory using default CNs, since very large errors between predictions and observations were achieved. Therefore, the SCS-CN model was adjusted by manual trials tuning both λ and CN parameters until the maximum coefficient of efficiency E (see below) was achieved using optimal λ and CN (Table 3).

**Table 3.** Optimal values of the Soil Conservation Service-curve number (SCS-CN) model parameters used for runoff predictions in plots subject to prescribed fire and wildfire (Lezuza and Liétor, Castilla La Mancha, Spain).


Note: \* indicates the CN value of the first modeled event.

#### *2.5. Evaluation of Model Prediction Accuracy*

The runoff simulations of the linear regression equation and SCS-CN model were analyzed for "goodness-of-fit" with the corresponding observations. More specifically, the observed and simulated runoff volumes were visually compared in scatterplots. Then, the main statistics and the indexes of goodness-of-fit commonly used in literature (e.g., [27,48–50]) were adopted (Table 4): (i) the maximum, minimum, mean and standard deviation of both the observed and simulated values; (ii) the coefficient of determination (r2); (iii) the coefficient of efficiency of Nash and Sutcliffe [51] (E); (iv) the root mean square error (RMSE) and (v) the coefficient of residual mass (CRM, also known as "percent bias", PBIAS). Table 4 reports the equations and the range of variability of these indexes [52–55]. Generally speaking, these indexes are based on the analysis of the errors between simulations and predictions of the modeled hydrological variables.

**Table 4.** Indexes and related equations and range of variability to evaluate the runoff prediction capacity of the linear regression and curve number models in forest plots subject to prescribed fire and wildfire (Lezuza and Liétor, Castilla La Mancha, Spain).


Notes: *n* = number of observations; Oi, Pi = observed and predicted values at the time step i; O = mean of observed values.

#### **3. Results and Discussion**

#### *3.1. Hydrological Characterization*

#### 3.1.1. Wildfire-Affected Forest (Liétor)

During the observation period, only nine events (total rainfall of 413 mm) produced surface runoff. For these events, precipitation height and mean intensity were in the range 11.6–93.8 mm and 0.98–28.0 mm/h, respectively. Expectedly, in the burned plots the runoff (on average 0.60 mm with a maximum value of 2.20 mm) was higher compared to the unburned soils (average of 0.03 mm and maximum of 0.08 mm). This may be due to the reduced infiltration and some combination of sealing, soil water repellency, loss of surface cover and decrease in soil aggregate stability, for the loss of organic matter [61].

In the mulched soil the mean and maximum runoff was 0.53 and 1.65 mm, respectively (Figure 2). These volumes were lower compared to the runoff generated in the burned plots. This shows the effectiveness of mulching as post-fire management technique to reduce the runoff generation capacity of the burned soils. These results confirm several other studies about the efficacy of mulch application, in order to control the hydrological response of soil after wildfires (e.g., [8,62–66]).

**Figure 2.** Rainfall and runoff volumes observed in forest plots subject to prescribed fire (Lezuza, Castilla La Mancha, Spain). Vertical bars are the standard deviations.

A sudden increase in the runoff generation capacity was evident in the first event after the fire (21 October 2016), presumably due to the ash release (that sealed the soil) and changes in the physicochemical properties (as the depletion in the organic matter content, which reduces the aggregate stability of the soil) [28,67]. A temporal reduction in runoff generation was found in burned soils (treated or not). This indicates a decrease in the hydrological response over time since the fire, also noticed by several authors in the early storms immediately after wildfire (e.g., [68–70]). The higher runoff is due to both the changes in soil hydrological properties and to the reduction of the vegetal cover after fire. As a matter of fact, the development of a water-repellent layer (also due to the ash released by fire) over the soil surface and the destruction of soil aggregates reduce water infiltration and thus increase runoff [71,72]. Over time, the shrub and herb vegetation quickly recovery, which decreases the runoff generation on the soil left bare by wildfire [73].

#### 3.1.2. Forest Subjected to Prescribed Fire (Lezuza)

Sixteen storms (totaling 368 mm) produced runoff. The mean runoff from these rainfall events was 0.39 mm and the maximum was 0.69 mm. In the burned plots, the mean and maximum runoff volumes recorded were 0.40 and 0.75 mm, respectively (Figure 3). For few events, runoff from burned soils was lower compared to the control plots, while, for the majority of the monitored storms, the soil subjected to prescribed fire produced noticeably more runoff compared to the unburned plots. This waiving soil response to storms confirms the low impacts of low-intensity fires on the hydrological response of soils, already observed by several authors (e.g., [28,74]). This means that the prescribed fire has a limited potential to change the soil properties that drives the hydrological behavior, such as the repellency and infiltrability. However, as observed for the wildfire, attention should be paid to the first rainfall events occurring immediately after fire, when the removal of almost all the vegetal cover (including the litter) may leave the soil bare and thus exposed to rainfall erosivity.

**Figure 3.** Rainfall and runoff volumes observed in forest plots subject to wildfire (Liétor, Castilla La Mancha, Spain). Vertical bars are the standard deviations.

#### *3.2. Hydrological Modeling*

#### 3.2.1. Linear Regression

The simulation of runoff volumes gave satisfactory results for the unburned soil both in Lezuza and Liétor, as shown by the values of E (close or over 0.75) and r<sup>2</sup> (>0.62) indexes and the closeness between predictions and observations (mean error of 10-15%). RMSE values (0.01, Lezuza, and 0.07, Liétor) were under the limit of acceptance of Table 4 (half std. dev.) only for the runoff measured in Lezuza (0.165), but not in Liétor (0.005). In general, the linear models tended to overestimate the runoff volumes in unburned and burned and mulched plots (CRM < 0), while underestimating the runoff in burned plots (CRM > 0; Figure 4). In the latter condition, the maximum runoff values were noticeably underestimated (difference between 25 and 40%; Table 5).

The performance of the linear regression Equation (7) was only acceptable but not satisfactory in burned soils (with or without treatment), because E (between 0.52, soil burned by wildfire, and 0.62, soil burned by wildfire and then mulched) was lower than the suggested limit of Table 4 and the differences between the maximum values of observations and predictions were over 20%; only the mean values were close each others (error <12%; Table 5 and Figure 5a,b); the values of RMSE, which was acceptable in Lezuza (0.12 vs. a limit of 0.15), were 0.62 (burned and untreated plots) and 0.46 (burned and mulched soils) and therefore were over the acceptance limit (Table 4). This limited performance is mainly due to inaccurate prediction of the most intense rainfall–runoff event (21 October 2016), immediately following the fire. Moreover, for the burned soils, the RMSE values were higher than 50% of the standard deviation of observed runoff and thus not satisfactory; for soils burned by wildfire, also the coefficient of determination was poor (r<sup>2</sup> < 0.39). This unsatisfactory model performance is also visually shown by the large scattering of the simulations around the regression line (Figure 5a,b), which highlights a particular prediction inaccuracy for the first rainfall event (21 October 2016 recorded

in Liétor). This inaccuracy is due to the soil changes induced by the wild fire (e.g., water repellency, decreases in infiltration and interception), which alters its hydrological response to rainstorms, but disappear some weeks or a few months after the fire [7–9,26]. Moreover, the linear equation was not able to simulate the variability of the hydrological processes with the precipitation, since the same runoff was observed for the same precipitation. This means that linear regressions are not able to simulate with reliability the surface runoff produced in burned plots, although these models may give an indication at least of the magnitude of the hydrological response of soils under different precipitation input and conditions.

&RQWURO %XUQHG %XUQHGDQGPXOFKHG

**Figure 4.** Linear regressions between observed rainfall and runoff in plots subject to prescribed fire and wildfire (Lezuza, **Left**, and Liétor, **Right**; Castilla La Mancha, Spain). V is the runoff volume and P is the rainfall height, while r2 and p are the coefficient of determination and the significance level, respectively.

**Table 5.** Statistics and indexes to evaluate the runoff prediction capacity of linear regression models in forest plots subject to prescribed fire and wildfire (Lezuza and Liétor, Castilla La Mancha, Spain).


Notes: r2 = coefficient of determination; E = coefficient of efficiency; CRM = coefficient of residual mass; RMSE = root mean square error.

**Figure 5.** Scatterplots of runoff observations vs. predictions using linear regressions in plots subject to prescribed fire (Lezuza, **a**) and wildfire (Liétor, control soils, **b**, and burned soils **c**) (Castilla La Mancha, Spain).

#### 3.2.2. SCS-CN Model

The predictions of runoff volume became more accurate for the majority of fire (in terms of severity) and soil conditions (Figure 6a,b). Compared to the linear regressions, the runoff predictions improved in the unburned plots of Lezuza (E = 0.87 and r<sup>2</sup> = 0.92), but slightly worsened in all the plots of Liétor (E = 0.88 and r<sup>2</sup> = 0.95; Table 6).

(b) (c)

**2EVHUYHGUXQRIIPP**

**Figure 6.** Scatterplots of runoff observations vs. predictions using the curve number model in plots subject to prescribed fire (Lezuza, **a**), and wildfire (Liétor, control soils, **b**, and burned soils, **c**) (Castilla La Mancha, Spain).

Conversely, the runoff was predicted with greater accuracy in burned soils (with or without treatment), as shown by E > 0.80—with peaks of 0.96-0.97 after wildfire—and r2 > 0.94, these indicators being noticeably over the acceptance limit (Table 4). The RMSE values were always lower than 50% of the observed standard deviation (Table 6). In general, the SCS-CN model always showed a runoff overestimation (see CRM < 0). It should be noticed that the mean runoff values were predicted

with less accuracy compared to the linear regression equations (differences between predictions and observations over 20% in some cases), but the SCS-CN model was much more reliable in predicting the maximum runoff volumes (differences lower than 7% for soils burned by wildfire and 15% for plots subjected to prescribed fire; Table 6). This indicates that, when accurate predictions of runoff are required to control the flooding risk (linked to the highest runoff volumes), the SCS-CN should be preferred to the simpler linear regressions.

**Table 6.** Statistics and indexes to evaluate the runoff prediction capacity of the curve number model in forest plots subject to prescribed fire and wildfire (Lezuza and Liétor, Castilla La Mancha, Spain).


Notes: r2 = coefficient of determination; E = coefficient of efficiency; CRM = coefficient of residual mass; RMSE = root mean square error.

Some additional considerations about SCS-CN model application in the experimental conditions should be made.

First, very low values of CN and λ were provided in this study as input to the model, in order to predict with accuracy runoff after the two fire-severity conditions. This means that the water losses during and immediately after the rainfall (reflected by *Ia* and S, such as water storage in the soil dips, interception, infiltration and evapotranspiration) are very high and the storms produce very small runoff volumes. In more detail, the small CN simulates a large water storage capacity (S) of soil through the infiltration process and λ must be decreased even by three orders of magnitude to simulate the very low initial water losses, due to interception and evapo-transpiration.

Second, both for wildfire and prescribed fire, unrealistic input parameters are required to simulate such a minimal runoff generation capacity of these soils. As a matter of fact, values of 15–16 (after wildfire) or even 0.25–3 (after prescribed fire) for CN and 0.0001 for λ against common values over 30 for CN and 0.2 for λ are needed to fit the runoff predictions to the corresponding observations. This should be taken into account when the SCS-CN model must be implemented in soils having a small hydrological response.

Third, a unique CN value as input for the SCS-CN model is not able to reproduce the increase in runoff immediately after wildfire. The worsening of the hydrological response of the burned soil both after wildfire and prescribed fire has been shown by a number of studies (e.g., [9]) and particularly in Mediterranean forests (e.g., Keizer et al. [8], in a Portuguese eucalypt forest; Lucas-Borja et al. [7,62], in pine forests of Central Spain). This increase is mainly due to soil water repellency and vegetation cover removal due to fires, but these effects disappear some months after fire. In order to simulate the hydrological effects of a repellent and almost bare soil, it is necessary to increase the CNs in this so-called

"window-of-disturbance" [15,75] up to values that may be noticeably high. Soil mulching "smoothes" the increase in the soil hydrological response and this requires a lower increase in CN values.

#### **4. Conclusions**

This study evaluated the runoff prediction capacity of the SCS-CN model in comparison with linear regression equations after storms in two pine forests of Central-Eastern Spain affected by wildfire (with or without a rehabilitation treatment using straw mulching) or prescribed fire.

The simulation of runoff volumes by the linear regression gave satisfactory results only for the unburned soils. Conversely, for the burned plots, the linear regressions failed in simulating the runoff with reliability. The SCS-CN model was instead accurate to predict the runoff volume particularly in burned soils, also when mulching was applied. Although the mean runoff was predicted with less accuracy compared to the linear equations, the model performance was very satisfactory in predicting the maximum volumes. Moreover, all the soil conditions (unburned, burned and burned and mulched) were simulated with reliability. To reproduce the peculiar hydrology of the experimental areas, very low values of CNs and initial abstraction were required, which may appear unrealistic; moreover, the post-fire hydrological window-of-disturbance could be reproduced only by increasing the CN for the storms occurring few months after wildfire.

The performances of the two tested models indicate that, in Mediterranean forests subject to the fire risk, the use of simple linear equations is suggested for predicting runoff generated by relatively low storms, while the SCS-CN model is more reliable and therefore advisable when runoff predictions are needed to control the flooding risk.

Overall, the study has confirmed the viability of the SCS-CN method to reproduce the complex hydrological response of unburned and burned (and treated or not) soils of Mediterranean forests. Although this assumption is limited to the experimental conditions, the results are encouraging towards larger applications of this method in other climatic and geomorphologic conditions. However, further modeling studies are needed, in order to explore the runoff prediction capacity of the model in fire-affected forests with different ecological and soil characteristics. These studies should also be enlarged from the plot to the watershed scale, using more complex hydrological models based on the SCS-CN method. Once validated in a wide range of environmental contexts, the use of these models may support the land managers to control runoff and erosion in mountain forests that are prone to both the wildfire and hydrogeological risks.

**Author Contributions:** Conceptualization, D.A.Z.; methodology, A.L., D.D., P.D., P.A.P.-A. and B.G.C.; formal analysis, G.B. and M.E.L.-B.; investigation, A.L., D.D., P.D., P.A.P.-A. and B.G.C.; resources, M.E.L.-B.; data curation, A.L., D.D., P.D., P.A.P.-A. and B.G.C.; writing—Original draft preparation, D.A.Z. and M.E.L.-B.; writing—review and editing, D.A.Z., G.B. and M.E.L.-B.; supervision, D.A.Z., G.B. and M.E.L.-B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Pedro Antonio Plaza Álvarez was supported by a predoctoral fellowship from the Spanish Ministry of Education, Culture and Sport (FPU16/03296).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

# **Development of a Parametric Regional Multivariate Statistical Weather Generator for Risk Assessment Studies in Areas with Limited Data Availability**

#### **Saddam Q. Waheed 1,2,\*, Neil S. Grigg <sup>1</sup> and Jorge A. Ramirez 1,**†


Received: 10 July 2020; Accepted: 7 August 2020; Published: 11 August 2020

**Abstract:** Risk analysis of water resources systems can use statistical weather generators coupled with hydrologic models to examine scenarios of extreme events caused by climate change. These require multivariate, multi-site models that mimic the spatial, temporal, and cross correlations of observed data. This study developed a statistical weather generator to facilitate bottom-up approaches to assess the impact of climate change on water resources systems for cases of limited data. While existing weather generator models have impressive features, this study suggested a simple weather generator which is straightforward to implement and can employ any distribution function for variables such as precipitation or temperature. It is based on (1) a first-order, two-state Markov chain to simulate precipitation occurrences; (2) the use of Wilks' technique to produce correlated weather variables at multiple sites with the conservation of spatial, temporal, and cross correlations; (3) the capability to vary the statistical parameters of the weather variables. The model was applied to studies of the Diyala River basin in Iraq, which is a case with limited observed records. Results show that it exhibits high values (e.g., over 0.95) for the Nash–Sutcliffe and Kling–Gupta metric tests, preserves the statistical properties of the observed variables, and conserves the spatial, temporal, and cross correlations among the weather variables in the meteorological stations.

**Keywords:** statistical weather generator; stochastic process; Diyala River basin; Wilks' technique

#### **1. Introduction**

Climate change impacts are of increasing concern to hydrologists who assess risks in the management of water resources systems. Their models of climate scenarios for extreme events can be derived from global climate models (GCMs), stochastic-statistical weather generators (SWGs), or a combination. Although they have their own advantages, some argue that the GCM scenarios are inadequate and limit decision-making options because they represent only specific scenarios for climatic variability and have large uncertainties [1–6]. On the other hand, others think that SWGs can produce a wide range of scenarios to study system responses and provide more insights about the system performance under climate change [7–10]. The drawbacks of the SWGs are that they have a stochastic-basis and cannot provide future change insights. Therefore, the SWGs and GCMs have been linked to generate forecasting scenarios and to assign a probability of each SWG scenario by fitting a distribution to the GCM outcomes [11–13]. In this way, SWGs can then be used to generate probabilistic synthetic scenarios with the aid of the GCM information and which are statistically similar to observed data and used to investigate which climate states cause system failure [4,14–23]. Where historic records are limited, synthetic weather sequences based on SWGs are especially suitable [24].

Given the previous work, the main objective of this paper is to develop a SWG that can be used in a bottom-up approach to generate daily synthetic scenarios to evaluate the impacts of long-term climate change on system performance and suggest robust adaptations to cope with anticipated negative impacts that will be examined in a follow up study. Emphasis is placed on areas with low data availability, and the model is demonstrated for Diyala River basin in Iraq for the four historic weather variables (e.g., precipitation, maximum and minimum temperatures, and wind speed magnitude) with daily time steps from 1948 to 2006.

#### **2. Literature Review**

Generally, SWGs can be grouped into parametric, non-parametric, and semi-parametric methods. In the parametric method, the weather variables are assumed to fit one continuous probability distribution or two combined distributions. The parameters are usually estimated from historic observations [24–28]. In the non-parametric method, the weather variables are resampled from historic observations using techniques such as empirical distributions, neural networks, and maximum entropy bootstrap [29–32]. The semi-parametric method is a mixture between parametric and non-parametric methods.

Albeit other approaches have their advantages, the parametric SWG in the bottom-up approach is preferable because the parameters can be altered to simulate different weather scenarios and facilitate climate change studies [16]. Verdin et al., [15], Furrer and Katz [18], Buishand and Brandsma [33], Seneviratne et al., [34] noted that the non-parametric method has limitations in generating extreme events because values can only be in the range of the observations. Using only the observed sequences ignores climate change's impacts on altering the intensities of the variables and is insufficient in assessing the future response of water resources systems because it leads to single results corresponding only to these observed sequences [22,23,25,34–36].

Most existing SWGs are for single sites and cannot capture the spatial and cross correlations between the variables, which are essential for generating realistic climate change scenarios. Schaake et al., [37] stated that "relationships between physically dependent variables like precipitation and temperature should be respected". Single site SWGs can fail to capture the extreme events of the generated runoff, which are essential to develop realistic adaptation strategies to cope with flood and drought events, especially where a high runoff in one sub-basin can be offset by the low runoff in adjacent sub-basins [26,35,38].

Moreover, the misrepresentation of spatial and cross correlations (e.g., correlations between the precipitation and temperature) leads to biased generated streamflows as this correlation determines the water availability for evapotranspiration and snowmelt [32,39,40]. Therefore, SWGs should capture the characteristics of each site and the spatial dependence among them.

Recently, multi-site and multi-variable SWGs have been developed using different approaches. Steinschneider and Brown, [4] developed a semi-parametric model using a k-nearest-neighbor resampling scheme to simulate multiple spatially distributed variables using wavelet decomposition and autoregressive model to account for low-frequency oscillations. They used a Markov chain of first-order with three states to identify the precipitation states (e.g., dry, wet, and extremely wet). This model had difficulty in preserving the weather statistics besides the cross correlation. Additionally, it is not clear how to diagnose the differences between the precipitation states (e.g., wet and extremely wet).

Srivastav and Simonovic, [32] developed a non-parametric model using the maximum entropy bootstrap technique to capture the time-dependent structure and statistical characteristics. They used an orthogonal transformation to capture the spatial correlations. Even though the model preserves the historical characteristics, Verdin et al., [15] and Chen et al., [40] showed that the maximum entropy bootstrap technique is limited to the historical data range leading to inadequacy in climate change studies. It is difficult to employ this model to create different climate scenarios through variations of parameters.

Li and Babovic, [41] proposed a two-stage parametric model using an empirical copula to generate spatial distribution templates. Then, they developed a rank ordering technique that depended on historic data ranks with an empirical copula technique to preserve the correlations between the variables. The model preserves correlations between the variables and sites but is limited to the historic record length. For example, the model cannot generate more than 30 years of simulation if the historic observations are 30 years. Therefore, the model is not useful in areas with limited data length as an insufficient projection length may lead to wrong conclusions in risk assessment studies [42,43].

Verdin et al., [15] presented a model using a Bayesian hierarchical technique. The precipitation amounts are modeled using gamma distributions and maximum and minimum temperatures are modeled using a normal distribution. The statistical coefficients within them are modeled as spatial Gaussian processes to account for the correlations. Besides the complexity of model structure, the model has difficulty in preserving the statistical properties of the variables (especially the standard deviation of the minimum temperature is extremely underestimated by the model). Additionally, the model underestimates the spatial correlation between the variables. Furthermore, their results do not demonstrate the model's ability to preserve the cross correlation between the variables as well as the temporal correlation.

#### **3. Model Description**

The goal here is to develop a parametric regional weather generator (PR-WG) to generate daily stochastic weather variables that preserve their statistical parameters, such as the mean and standard deviations, as well as the spatial, temporal, and cross correlations among them. It should be easy to implement and adapt by altering the statistical parameters to generate synthetic future climate scenarios. The generated scenario must exceed the historic record length and observation range.

The novel contribution is to use a parametric approach to create a flexible model that can adapt to any continuous probability distribution. This will enable the use of the most accurate distribution for each weather variable, and the user can employ other distributions according to the data availability and scope of the study.

#### *3.1. Precipitation States*

The first step in developing the PR-WG is to establish the precipitation states. They are defined here as: wet days if the daily amounts equal or exceed 0.1 mm and dry days otherwise. This is similar to the approach by Verdin et al., [15] and Li and V. Babovic [41]. The approach is to use the first-order two-state Markov chain (FTMC), which is the most popular method to produce dry and wet precipitation occurrences. It works well in different climate types and performs as well as higher Markov chain orders [21,22].

Let *S*(*k*,*t*,*m*) denote the precipitation state (*S* = 0 is a dry day and *S* = 1 is a wet day) at spatial location *<sup>k</sup>* <sup>∈</sup> <sup>N</sup>, time index *<sup>t</sup>* <sup>∈</sup> <sup>N</sup> in days, and month index *<sup>m</sup>* <sup>=</sup> {1,2, ... 12}. The dry or wet day occurrence is obtained from the following conditional probabilities:

$$\Pr\left(S\_{(k,t,m)}=0\mid S\_{k,t-1,m}=0\right) = \kappa\_0 \; ; \Pr\left(S\_{(k,t,m)}=1\mid S\_{k,t-1,m}=0\right) = 1-\kappa\_0\tag{1}$$

$$\Pr\left(S\_{(k,t,m)}=1\mid \mathbb{S}\_{k,t-1,m}=1\right) = \kappa\_1; \Pr\left(S\_{(k,t,m)}=0\mid \mathbb{S}\_{k,t-1,m}=1\right) = 1-\kappa\_1\tag{2}$$

where, κ<sup>0</sup> is the probability of a dry day following a dry day, and κ<sup>1</sup> is the probability of a wet day following a wet day. These probabilities were estimated from the daily historical precipitation observations for each month.

#### *3.2. Precipitation Amount*

Precipitation amounts were calculated by using the joint probability distribution between the occurrence and amount. For example, once a wet day is predicted from the FTMC, the precipitation amount is calculated. A skewed normal distribution (SN) was selected because it was recommended by other researchers and estimates the daily precipitation amount better than other distributions such as exponential, gamma, Weibull, mixed-exponential, and generalized Pareto in capturing the mean, standard deviation, and extreme values [20,21,36,44,45].

Let P denote the precipitation amount in mm/day and ঌ<sup>Ψ</sup> denote the indicator of precipitation state condition ψ. P returns to a value obtained implicitly from Equation (4) [46] if the condition ψ holds (ঌ[*S*=1]) and returns to zero otherwise ঌ[*S*=0] , as follows:

$$P\_{(k,t,m)} = \begin{cases} \text{SN} \left( \mu\_{P\prime} \right. \left. \sigma\_{P\prime} \left. \gamma\_P \right) \right. & \left. for \begin{array}{l} \mathbb{I}\_{[S(k,t,m)=1]} \\ for \begin{array}{l} \mathbb{I}\_{[S(k,t,m)=0]} \end{array} \end{array} \tag{3}$$

$$\theta\_{(k,t,m)} = \frac{6}{\mathcal{V}\_p(k,m)} \left\{ \left[ \frac{\mathcal{V}\_{p(k,m)}}{2} \left( \frac{P\_{(k,t,m)} - \mu\_{p(k,m)}}{\sigma\_{p(k,m)}} \right) + 1 \right]^{\frac{1}{5}} - 1 \right\} + \frac{\mathcal{V}\_{p(k,m)}}{6} \tag{4}$$

where θ is the matrix of the standard normal deviates θ ~ N(0,1) R, and μ*p*, σ*p*, and γ*p*, are the mean, standard deviation, and skew coefficient of the precipitation for month *m*. The values of the parameters μ*p*, σ*p*, and γ*<sup>p</sup>* were estimated from the daily historical observations using the method of maximum likelihood estimation (MLE).

#### *3.3. Maximum and Minimum Air Temperature*

The maximum and minimum daily air temperatures are usually modeled by the normal distribution (N) [47,48]. Let *TX* and *TN* denote the maximum and minimum daily air temperature in ◦C, respectively. In which, *TX* is (and *TN*) is:

$$T\_{X(k)} \sim N(\mu\_{X(k)\prime}, \sigma\_{X(k)}) \tag{5}$$

where μ*<sup>X</sup>* and σ*<sup>X</sup>* are the mean and standard deviation of TX, respectively. Solving Equation (5) for each month m according to ঌ<sup>Ψ</sup> (to account for precipitation state effects), Tx and TN can be computed as:

$$T\_{X(k,t,m)} = \mu\_{x0(k,m)} + \sigma\_{X0(k,m)} \times \mathbb{I}\_{(k,t,m)} \text{ for } \mathbb{I}\_{[S(k,t,m)=0]} \tag{6}$$

$$T\_{X(k,t,m)} = \mu\_{\mu x 1(k,m)} + \sigma\_{X1(k,m)} \times \mathbb{I}\_{(k,t,m)} \text{ for } \mathbb{I}\_{[S(k,t,m)-1]} \tag{7}$$

$$T\_{N(k,t,m)} = \mu\_{\mu N0(k,m)} + \sigma\_{N0(k,m)} \times \delta\_{(k,t,m)} \text{ for } \mathbb{I}\_{[S(k,t,m)=0]} \tag{8}$$

$$T\_{N(k,t,m)} = \mu\_{\mu N1(k,m)} + \sigma\_{N1(k,m)} \times \delta\_{(k,t,m)} \text{ for } \mathbb{I}\_{[S(k,t,m)=1]} \tag{9}$$

where, μ*X*0, μ*X*1, μ*N*0, μ*N*1, σ*X*0, σ*X*1, σ*N*0, and σ*N*<sup>1</sup> are the monthly mean and standard deviation for the maximum and minimum air temperature (◦C/day) for *S* = 0 and 1, respectively, and υ and δ are the matrices of standard normal deviates, such that υ and δ ~ N (0,1) R. The parameter values of Equations (6)–(9) were estimated from the historic observations using MLEs.

#### *3.4. Wind Speed Magnitude*

Ref. [49] showed that the most accurate function to simulate the daily wind speed magnitude (WS) is Weibull with three and two parameters, respectively, followed by gamma. Given the condition that wind speed is affected by precipitation states and amounts [50], the selected distribution must be decomposed into the same distribution type. As the Weibull distribution cannot be decomposed into two Weibulls (although gamma can be [51]), wind speed magnitude was modeled by the gamma distribution (GM) in this study. Let WS denote the daily wind speed magnitude (m/s) for k locations, as follows:

$$\text{WS}\_{(k)} \sim \text{GM}\left(\alpha\_{(k)}, \beta\_{(k)}\right) \tag{10}$$

where α and β are the shape and scale parameters, respectively. Similarly for the temperature, the WS for each month m, according to ঌΨ, was estimated implicitly from the following equations:

$$\Lambda\_{\left(k,t,m\right)} = \frac{\beta\_{0\left(k,m\right)} - a\_{0\left(k,m\right)}}{\Gamma\left(a\_{0\left(k,m\right)}\right)} \int\_{0}^{\text{WS}\_{\left(k,m\right)}} h^{a\_{0\left(k,m\right)}-1} \exp^{-h/\beta\_{0\left(k,m\right)}} \, d\mathbf{l} \, f \, \text{or} \, \mathbb{I}\_{\left[\mathbf{S}\left(k,t,m\right)=0\right]}\tag{11}$$

$$\lambda\_{(k,t,m)} = \frac{\beta\_{1(k,m)} ^{-\alpha\_{0(k,m)}}}{\Gamma\left(\alpha\_{1(k,m)}\right)} \int\_0^{\text{WS}\_{(k,m)}} h^{\alpha\_{1(k,m)}-1} \exp^{-h/\beta\_{1(k,m)}} \, dl \, f \, \text{or} \, \mathbb{I}\_{[\text{S}(k,t,m)=1]} \tag{12}$$

where α0, α1, β0, and β<sup>1</sup> are the shape and scale parameters for *S* = 0 and 1, respectively, for each month m, *h* is an independent parameter, and λ is the cumulative probability, which is distributed uniformly—λ~ U [0, 1], R. The shape and scale parameters were estimated from the historic observations using MLEs.

#### **4. Model Implementation**

The parametric SWG should conserve the spatial, temporal, and cross correlations of the historic observations of the four weather variables. The concept is to study the behavior of the variates θ, υ, δ, and λ, hereafter referred to as anomalies. The correlations between those anomalies should be identified so the generated weather values are statistically similar to the observed values and conserve spatial, temporal, and cross correlations. The implementation of the PR-WG consists of two stages, namely preprocessing and postprocessing, as shown in Figure 1.

#### *4.1. Preprocessing: Parameter Estimation and Matrix Preparation*

In order to specify the wet and dry occurrences, a random uniform variate y ~U(0, 1) must be drawn and compared with the transition probabilities obtained from Equations (1) and (2). For multi-site precipitation, the anomalies (referred to as Y R) that identify the states in k locations must be correlated so that the generated states S are correlated to the historic observations. Wilks' method was selected to generate correlated anomalies Y~ N(0,1) at multiple sites. It is simple and more efficient than hidden Markov and k-nearest neighbor methods [52], accurate in generating the correlations of monthly interstations [53], and the most cited method compared to other approaches [54].

Assume S (1,m) and S (2,m) are the precipitation states on month m at sites k = 1 and k = 2. To generate realistic sequences of the precipitation states at these two sites, the correlation (ω) between their corresponding anomalies Y, ω(1,2) = corr (Y(1,m), Y(2,m)) must be computed. The parameter ω was determined by generating different sets of Ý at the two sites with different arbitrary correlation values {ω´ 1, ω´ 2, ... }, ω´ <sup>1</sup> =corr (Ý (1,m), Ý (2,m)), identifying the precipitation states at the two locations S´ <sup>1</sup> and S´ 2, and calculating the corresponding correlation {ε´1, ε´2, ... }, ε´1(1,2) = corr (S´(1,m), S´ (2,m)). Then, a regression line between ε´ and ω´ sets was fitted to identify the relationship between them. Using this regression equation with the observed precipitation state correlation ξ, the parameter ω can then be found. A synthetic example is shown in Figure 2a, in which selecting a 0.858 correlation between the pair anomalies (ω) will produce 0.785 correlation between the pair states (ξ) at the two locations.

The process should be repeated for each station pair and lead to the number of realizations of k (k-1)/2 and be repeated for each month m to create the anomalies matrix <sup>ω</sup><sup>s</sup> <sup>∈</sup> <sup>R</sup>. The <sup>ω</sup><sup>s</sup> matrix is then used to develop Y that produces correlated precipitation states in k locations for month m, using the multivariate normal distribution as follows,

$$Y = f\left(\mu\_{y\_Y} \, \Sigma\right) = \frac{1}{\sqrt{\Sigma \left(2\pi\right)^d}} \exp\left(-\frac{1}{2}\left(y - \mu\_y\right) \Sigma^{-1}\left(y - \mu\_y\right)\right) \tag{13}$$

**Figure 2.** (**a**) An example of Wilks' technique for precipitation states; (**b**) and (**c**) are examples of Wilks' technique to obtain ϕ<sup>z</sup> for Tx and WS, respectively, for station k of month m.

The variable μ*<sup>y</sup>* denotes the 1-D mean vector for the anomalies *Y*, Σ denotes the covariance matrix, and d is an independent parameter. In this case, μ = [0, 0, ... , 0]k×<sup>1</sup> and the variance is 1, so the covariance matrix Σ<sup>s</sup> becomes the correlation matrix ωs.

The matrix ω<sup>s</sup> must be a positive-definite matrix (e.g., the matrix is symmetric and all its eigenvalues are positive) to be implemented in Equation (13). Since the elements of ω<sup>s</sup> were calculated empirically, ω<sup>s</sup> is usually a non-positive matrix. Comparing to the work of others, the most precise method to obtain a positive-definite matrix is the iterative spectral with Dykstra's correction (ISDC) [55], as follows:


After generating the matrix S at k and m, the next step is to simulate the weather variables (e.g., P, TX, TN and WS). The idea here is to examine the anomalies of these variables and generate the weather variables with the same observation properties. To account for all the spatial and cross correlation between the variables, their anomalies (θ, υ, δ, and λ) must be correlated. The temporal correlation, identified by the Lag-1 day auto-correlations, for the TX, TN and WS must also be considered. Since the precipitation amount is an intermittent variable, the auto-correlation is not considered. The following procedure was suggested to achieve this purpose. First, arrange the weather variable matrix V as follows,

$$
\begin{bmatrix}
V\_{1,1}^1 & V\_{1,2}^1 & \dots & V\_{1,k}^n \\
V\_{2,1}^1 & V\_{2,2}^1 & \dots & V\_{2,k}^n \\
\vdots & \vdots & \dots & \vdots \\
V\_{t,1}^1 & V\_{t,2}^1 & \dots & V\_{t,k}^n \\
\end{bmatrix} \tag{14}
$$

where, *V* represents the observed weather variable value and n denotes the weather variable rank (P, TX, TN and WS), *n* = {1, 2, 3, 4}. The total number of the rows will be T = month days × year numbers, the columns will be K × N, and the aisle will be M. This matrix arrangement enables us to consider all the spatial and cross correlations between the weather variables. Next, extract the anomalies matrix Z <sup>∈</sup> <sup>R</sup> from V using Equations (3) and (4) for P; Equations (6)–(9) for Tx and TN; Equations (11) and (12) for the WS, after estimating their parameters (e.g., μp, σp, γ<sup>p</sup> for P, μX0, μX1, σX0, σX1 for Tx, μN0, μN1, σN0, σN1 for TN, and α0, α<sup>1</sup> β0, β<sup>1</sup> for the WS).

The Z matrix represents the anomalies of the weather variables and their elements have spatial, cross-, and auto-correlation magnitudes. To generate the Z matrix with the same observation properties, these correlations must be preserved. The first step done here was to estimate autoregressive model of order 1, AR(1), coefficients for the anomalies (ϕz) so that the generated variables have the observed AR(1) value (ϕv) applying Wilks' technique. For illustration, synthetically assume that the values of μX0, μX1, σX0, σX1 are 11.72, 9.12, 3.71, 2.21 (Co/day), respectively, and ϕ<sup>v</sup> is 0.82 at station k of month m. The adopted procedure for obtaining the ϕ<sup>z</sup> is as follows:

1) Generate the standard normal random deviate set y;y~N (0,1).

⎡

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣


This procedure must be done for all Tx and TN of each k and m. For the WS, the procedure is the same except for Step 5, converting x so it is uniformly distributed to get the WS anomalies. For example, let us assume that α0, α1, β0, and β<sup>1</sup> are 4.04, 3.22, 0.62, 0.71, respectively, and the ϕ<sup>v</sup> is 0.54. The corresponding ϕ<sup>z</sup> will be 0.56, as shown in Figure 2c. This procedure allows us to preserve the auto-correlation of Tx, TN, and the WS.

The final step of the preprocessing stage is to construct the positive-definite correlation matrix of the variable anomalies ωV, as done for precipitation states using ISDC. Building the ω<sup>V</sup> allows us to preserve all the spatial, temporal, and cross correlations between the variables.

#### *4.2. Postprocessing Stage: Variable Generation*

After building all matrices and estimating the parameters in the preprocessing stage, the four weather variables can be generated for any time length of interest, as follows:


$$Z\_{\rm std(k)} = \frac{Z\_{(k)} - \mu\left(Z\_{(k)}\right)}{\sigma\left(Z\_{(k)}\right)}\tag{15}$$

where *Zstd* represents the standardized anomalies Z of Step 5, and μ(*Z*) and σ(*Z*) are the mean and standard deviation of *Z*, respectively.


$$Z\_{lI(k)} = 0.5 \times err\left(\frac{Z\_{(k)} - \mu\left(Z\_{(k)}\right)}{\sqrt{2}\ \sigma\left(Z\_{(k)}\right)}\right) + 0.5\tag{16}$$


#### **5. Case Study and Data**

The developed PR-WG was tested in the Diyala River basin, which is a transboundary basin between Iran and Iraq with a total stream length of 217 km and basin area of 16,760 km<sup>2</sup> above Derbendikhan Dam, as shown in Figure 3. In previous work, Waheed et al., [5] implemented the daily weather data (e.g., precipitation, maximum and minimum temperature, and wind speed) in this basin at a 0.5◦ spatial resolution from 1948 to 2006 and explained the implementation procedure. In this follow up study, the historic forcing data were used to validate the proposed PR-WG and test its performance. The reader should refer to the original paper for more details about the data implementation and their validation in the basin.

**Figure 3.** Diyala River basin in Iraq with grid-cell numbers.

#### **6. Results and Discussion**

#### *6.1. Model Performance Evaluation*

The PR-WG was tested for its daily performance with historic observations for the period between 1948 and 2006, e.g., 58 years, in a grid composed of 24 grid-cells. The Nash–Sutcliffe coefficient efficiency (NSCE; [56]) and the Kling–Gupta efficiency (KGE; [5]) were used to evaluate the PR-WG's ability to produce spatially correlated precipitation states S similar to the observed values, as follows:

$$\text{NSCE} = 1 - \frac{\sum \left( \text{Sim}\_{i} - \text{Obs}\_{i} \right)^{2}}{\sum \left( \mu\_{\text{sim}} - \text{Sim}\_{i} \right)^{2}} \tag{17}$$

$$\text{KGE} = 1 - \sqrt{\left(\frac{\mu\_{\text{sim}}}{\mu\_{\text{obs}}} - 1\right)^2 + \left(\frac{\sigma\_{\text{sim}}}{\sigma\_{\text{obs}}} - 1\right)^2 + \left(\rho - 1\right)^2} \tag{18}$$

where *Sim* and *Obs* are the simulations (e.g., the PR-WG outcomes) and the observations of the time index *t*, respectively; μ*obs*, σ*obs*, μ*sim*, and σ*sim* are the mean and standard deviation of the observations and simulations (e.g., the PR-WG outcomes), respectively, and ρ is the correlation coefficient between the observations and simulations.

Figure 4 shows the comparison of 10 separate daily simulations each of the same observation length (e.g., 58 years) of PR-WG monthly dry and wet occurrences in gray color dots. The average of these 10 simulations is calculated and plotted in blue dots. The A 1-1 line is also plotted to ease the comparison. It is evident that the model works well to produce the number of dry and wet days, with KGE and NSCE values of 0.97. This result demonstrates the ability of the FTMC to produce the precipitation states well [21,22]. Figure 5 shows a comparison of pairwise correlations of the daily precipitation states calculated for each calendar month. It can be seen that the correlations are captured well by the PR-WG. The overall KGE and NSCE values are 0.98 and 0.99, respectively.

**Figure 4.** Comparison of the daily precipitation states between the observations and simulations for all months and grid-cells.

**Figure 5.** Comparison of the daily precipitation state correlation between the observations and simulations for each month for all grid-cells.

Figure 6 demonstrates the PR-WG performance to produce the statistical parameters (e.g., mean, standard deviation and skewness) of the four weather variables. The comparisons were done on a daily basis at each month for the 24 grid-cells. A daily time step series of 1000 years was generated to reduce the sampling bias and uncertainty in the simulations. However, the daily means of all variables and the standard deviations for Tx, TN and WS were perfectly produced by the model (KGE ≈1), while σp, and γ<sup>p</sup> are reasonably preserved (KGE = 0.96 and 0.86; NSCE = 0.98 and 0.93). The slight discrepancies are due to the stochastic nature of the process [57].

*Climate* **2020**, *8*, 93

**Figure 6.** Comparisons of the daily statistic parameters of the observations and simulations. (**a**–**c**) are the mean, standard deviation, and skewness of P. (**d**–**g**) are the mean and standard deviation of Tx and TN. (**h**,**i**) are the mean and standard deviation of the WS.

Figure 7 shows the daily median values with 0.05 and 0.95 quantiles of the daily values in the bounded areas, and the inverse cumulative distribution function (CDF−1) of the observed and simulated weather variables for grid-cell number 9, which is located in the basin heart (see Figure 3). It is seen that PR-WG well preserves the daily medians for all months. Moreover, the quantile daily estimates show good agreement with the observation quantiles, proving the model's ability to capture the maximum and minimum daily weather values. It is also noticeable from the inverse CDF the observed and simulated weather values are very close, evincing the validity of the selected distribution types. Furthermore, the simulated daily values of quantile 1 exceed the observation values which demonstrates the model's ability to produce values beyond the observation ranges.

**Figure 7.** Comparisons of the daily observed and simulated values for the medians with daily 0.05 and 0.95 quantiles in the bounded areas, and the CDF−<sup>1</sup> for the four weather variables.

Figure 8 shows the spatial and cross correlation coefficient matrices of the observations and simulations for one month (e.g., m=1), while Figure 9 shows the spatial and cross correlation comparison for all variables for each m calculated at daily time steps. The number of columns of the V matrix (see Section 4.1) are 4 × 24 = 96. Therefore, the V dimensions are 96 × 96, in which the values are from 1 to 24 for P, 25 to 48 for Tx, 49 to 72 for TN, and 73 to 96 for the WS. It can be observed from Figure 8 that the observed correlation among the variables varies greatly across them. P and the WS are slightly less spatially correlated as compared with Tx and TN. These facts are in line with Srivastav and Simonovic [32] and Verdin et al. [15]. It is also noticeable from Figures 8 and 9 that the model preserves the spatial and cross correlation well among the variables. The overall KGE and NSCE values are 0.96 and 0.97, respectively.

**Figure 8.** Spatial and cross correlation coefficients of the daily observed (**a**) and simulated variables (**b**).

**Figure 9.** Spatial and cross correlation comparison of the daily weather variables for each month.

Figure 10 demonstrates the PR-WG capability to preserve the Lag-1 day auto-correlations of Tx, TN and the WS. It is noticeable that the values differ from month to month, they are less for the WS comparing to Tx and TN. However, the PR-WG captures these monthly variations very well regardless of their magnitudes with the overall KGE and NSCE values of 0.97 and 0.98, respectively.

The results presented here glimpse the model capability to preserve the statistical properties of the observations to synthesize the future scenarios. The proposed model demonstrated the Wilks technique ability to generate anomalies similarly to the observations. It is also seen that the hybrid structure of the AR and Wilks technique leads to generate data that preserve the temporal correlation beside the spatial and cross correlations.

*Climate* **2020**, *8*, 93

**Figure 10.** Lag-1 day auto-correlations of the weather variables Tx, TN, and WS, respectively, for all months.

The key advantage of PR-WG is that it is built to be a general model through studying the observation anomalies and mimicking them. Therefore, the model is anticipated to work well in different climate zones and topographies regardless of the data spatial and temporal scale. The model framework is flexible enough for locations observe short-term and long-term variations. Moreover, the user can reduce the cycle data length to meet their scope. e.g., they can use a data window of two weeks (or a week) instead of the monthly window that was used in this study. The computational expensive of implemented the pre-processing stage has to carefully examined.

#### *6.2. Model Validation*

In some cases, the proposed SWG produces negative daily values for precipitation. [58] indicated that the SN is not suitable when the skewness is greater than 4.5. However, in the study area, values of the skewness have not exceeded 4.5 (see Figure 6c), therefore the SN is applicable. The negatives of the daily values were checked and found to be less than 3% of the whole 1000-year time series in the 24 grid-cells. The suggestion of [32] to round the negative values to zero was considered, but it would affect the number of wet and dry calculations and the statistical parameters of precipitation. Instead, the negatives were rounded to 0.1 mm/day, which is assumed to be the minimum precipitation amount (see Section 3.1). This correction approach for negative values illustrates the slight differences in the simulated σ<sup>p</sup> and γ<sup>p</sup> (see Figure 5b,c). The user could apply another distribution function in cases where the SN is not applicable such as mixed-exponential [59,60], log-normal, gamma, etc. The key advantage of PR-WG is its flexibility in adopting any distribution of interest, such as these.

The second validation was done by checking if TN is greater than Tx and was found to be less than 1% of the whole 1000-year time series in the 24 grid-cells. [41] suggested to force Tx to be greater than TN through setting TN equal to Tx minus 1. This procedure will affect the auto-correlation of the TN. Instead, the Chen et al., [57] approach was applied as follows, if Tx < TN,

$$T\_{\mathcal{X}(k,t,m)} = T\_{\mathcal{N}(k,t,m)} + (\mu\_{\mu\mathbf{x}(k,m)} - \mu\_{\mu\mathbf{N}(k,m)}) + \sqrt{\sigma\_{\mathcal{X}(k,m)}^2 - \sigma\_{\mathcal{N}(k,m)}^2} \times z\_{\text{std}(k,t,m)} \text{ for } \mathbb{I}\_{[a\mathcal{X}\_{(k,t,m)} \ge a\mathcal{N}\_{(k,t,m)}]} \tag{19}$$

$$T\_{\mathcal{X}(k,t,m)} = T\_{\mathcal{N}(k,t,m)} + (\mu\_{\mu\mathcal{X}(k,m)} - \mu\_{\mu\mathcal{N}(k,m)}) + \sqrt{\sigma\_{\mathcal{X}(k,m)}^2 - \sigma\_{\mathcal{N}(k,m)}^2} \times z\_{\text{std}(k,m)} for \quad \mathbb{I}\_{[a\mathcal{X}\_{(k,t,m)} < a\mathcal{N}\_{(k,t,m)}]} \tag{20}$$

Equations (19) and (20) are conditioned on the precipitation states. For example, the σ and μ will turn to condition 0 if *S* = 0. In this case, the TX is guaranteed to be greater than TN and the auto, spatial, and cross correlations are preserved since they are multiplied by the anomalies *zstd*.

#### *6.3. Model Comparison*

For comparison purposes, two SWG models were selected to compare their performances with the PR-WG to further demonstrate the model applicability. The first model is the single site weather generator (WG) developed by Chen et al., [57] and Chen et al., [61]. The second is the two stages weather generator using an empirical copula (EC) approach developed by Li and Babovic [41]. To highlight the unique contribution of the PR-WG model, we focused on the model performances to maintain the spatial, cross, and temporal correlations. Figure 11 shows the daily performances of the WG and EC to for the spatial and cross correlations, while Figure 12 shows the temporal correlations for the sites and month. It is seen that the EC model works well in preserving the spatial, cross, and temporal correlation; the PR-WG is slightly superior to it. However, the KGE and NSCE for the spatial and cross correlations are 0.92 and 0.93, and for the temporal 0.95, and 0.96. It also notable that the WG model poorly preserves the spatial and cross correlations but has good ability to preserve the temporal correlation. This is because the model accounts for the temporal correlation only, where the simulated data were generated independently for all variables and sites which leads to poor spatial and cross correlation accuracy. The KGE and NSCE for the spatial and cross correlations are −0.29 and −8.3, and for the temporal 0.88, and 0.89. Although the EC approach works well in general, the only drawback is that its simulation time period must be identical to the historic observation, which prevents its usage in areas with limited data availability. This is because the post processing stage of the EC model employs a re-ranking technique that extracts the ranked variables directly from the historic observations. Therefore, the model length can only be the same as the historic observations, leading to less flexibility for future scenarios, especially in data scarce regions. The PR-WG has the advantage of producing the simulation length of interest, making it useful in areas with limited data availability besides maintaining the statistical characteristics.

**Figure 11.** Performance evaluation for empirical copula (EC) and weather generator (WG) models for preserving the spatial and cross correlations of the weather variables for each month.

**Figure 12.** Monthly performance evaluation for EC and WG models for preserving the temporal correlation (Lag-1 day) of the weather variables.

#### *6.4. Simulation of the Future Forecasting Scenarios*

The goal of the PR-WG is to be used later for climate variation assessments. The advantage of the model, besides the ability to preserve the statistical characteristics, is its flexibility to alter them to produce a wide range of different scenarios. However, defining the future scenario ranges to test a water resources system's performance in terms of the climate stress is a difficult task and dependent on many factors, including expert opinions [4].

These future scenarios will be applied in the Diyala River basin to discover the vulnerability of the Derbendikhan Dam and its reservoir. Moreover, different adaptation strategies will be suggested in order to test their capabilities to improve the system performance. Since the model is implemented on a stochastic basis, the future trend insights will be obtained from analyzing the GCM models. Then, this can be fed into the PR-WR to mimic the future trend as well as the statistical properties. For instance, multiplicative factors for the precipitation mean will be applied starting from a 0% change in the historical precipitation and annual linearly increasing (or decreasing) up to the specified value in the final period (e.g., +30% of the historical value). Forms other than the linear change can also be applied to synthesize the future forecast data.

#### **7. Conclusions**

It was shown that a PR-WG accurately preserves the statistical properties (mean, standard deviation, and skewness coefficient) of the weather variables (overall KGE and NSCE test values were 0.98). The PR-WG also preserves the spatial, temporal, and cross correlations among the weather variables. While other SWGs may have more features, the one developed in this study enables a bottom-up vulnerability assessment study to be implemented in areas with limited data availability.

The PR-WG effectively estimates the dry and wet day occurrences using a FTMC with overall KGE and NSCE values of 0.97, a result that is in line with those in [21,22]. The results also demonstrate the effectiveness of Wilks' technique to produce spatially correlated precipitation states (KGE of 0.98; NSCE of 0.99) and spatially and cross correlated weather variables (KGE of 0.96; NSCE of 0.97), as well as temporally correlated variables (KGE of 0.97; NSCE of 0.98). The model is also capable of preserving the maximum and minimum daily weather values as well as producing values beyond the observed ranges. Furthermore, the PR-WG outperforms the EC and WG models in preserving the spatial, cross, and temporal correlations in the meteorological stations.

While the PR-WG was validated in the Diyala River basin, it should be effective and applicable in other places and with other weather variables, such as solar radiation. The advantages of PR-WG are its flexibility to select any distribution function for each weather variable, ability to simulate any number of years within or beyond the historic observation length, capability to generate values outside the observation range, and its ability to produce synthetic scenarios through the alteration of the weather variable parameters for the study of climate change's impacts. The PR-WG is easy to construct and understand with little computational intensity to build the spatial and cross correlation matrices of the anomalies. Increasing computational power will facilitate the work.

**Author Contributions:** Conceptualization, S.Q.W., J.A.R., and N.S.G.; methodology, S.Q.W.; software, S.Q.W.; validation, S.Q.W.; formal analysis, S.Q.W., N.S.G., and J.A.R.; investigation, S.Q.W.; resources, S.Q.W.; data curation, S.Q.W.; writing—original draft preparation, S.Q.W.; writing—review and editing, S.Q.W., N.S.G., and J.A.R.; visualization, S.Q.W.; supervision, J.A.R. and N.S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Iraq Higher Committee for Education Development (HCED), grant number D1201077.

**Acknowledgments:** The authors are grateful to the Iraqi Ministry of Water Resources for assistance. The authors are also thankful to Xin Li and Vladan Babovic for their cooperation in implementing the EC model. The authors are also grateful to Colorado State University for providing its laboratories and supercomputer to run the model(s) and perform the analyses.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Hallegatte, S.; Shah, A.; Lempert, R.; Brown, C.; Gill, S. *Investment Decision Making under Deep Uncertainty-Application to Climate Change*; The World Bank: Washington, DC, USA, 2012.


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Climate* Editorial Office E-mail: climate@mdpi.com www.mdpi.com/journal/climate

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-5066-4