Assessment of Changes in Heatwave Aspects over Saudi Arabia during the Last Four Decades

: Heatwave (HW) number (HWN), frequency (HWF), duration (HWD), magnitude (HWM), and amplitude (HWA) are key aspects for interpreting and understanding HW characteristics world-wide. Most previous HW studies over the Kingdom of Saudi Arabia (KSA) focused only on the temperature extremes, so this study aims to assess the decadal changes, anomalies, and spatiotemporal variations in the ﬁve HW aspects over KSA during the last four decades (1982–2021) using the ClimPACT2 software. Daily gridded (0.25 ◦ × 0.25 ◦ ) maximum (TX) and minimum (TN) temperatures from the ECMWF-ERA5 reanalysis dataset were used to compute these heat wave (HW) aspects. The HW aspects were derived in ClimPACT2 using the Excess Heat Factor (EHF), the 90th percentile of TX (TX90), and the 90th percentile of TN (TN90), all based on the reference climate period of 1982–2011. The results showed that the decadal sum and anomaly of the ﬁve HW aspects increased gradually during the last four decades (1982–2021). The three indices showed that the maximum decadal sum of HWN (42 events), HWF (255 days), and HWD (145 days) occurred in the last decade. Additionally, the last decade has the maximum decadal sum of HWM (175–463 ◦ C) and HWA (189–471 ◦ C) as derived from TX90 and TN90, which is conﬁrmed by EHF, with ranges of 7–58 and 15–185 ◦ C 2 , respectively. Finally, the periods 2015–2021 and 1984–1986 recorded the highest and lowest values of annual HW aspects, respectively, across the study period.


Introduction
Climate change, global warming, and increasing trends in extreme weather events (heatwaves, droughts, flash floods, etc.) are significantly associated with human-caused factors such as greenhouse gases and aerosols [1,2].Global warming has increased steadily since the last decades of the twentieth century [3,4], and it is projected to continue increasing even under the lowest emission scenarios [5][6][7][8].This rise in global temperature leads to an increase in the different characteristics/aspects (number, frequency, duration, and intensity) of heatwaves (HWs) around the world [1,[9][10][11][12].The scientific research community agrees that an HW should be measured as a period of at least 2 or 3 consecutive days with high temperatures [13].HW is defined by the World Meteorological Organization (WMO) as a period of more than five consecutive days in which the daily maximum temperature exceeds its normal climate during the reference period (30 years) by 5 • C [14][15][16].HWs are usually specified based on the normal climate of the region.For example, a normal temperature in an area with a hot climate may be considered an HW in another region with a cooler climate.Moreover, the aforementioned criterion was not appropriate for varied climatic conditions and complex terrain (e.g., China) [17].This is why many climate researchers have used various definitions of HWs [18][19][20][21][22], and sometimes combine one or two factors such as maximum temperature, relative humidity, dew point temperature, and water vapor pressure [23][24][25][26][27][28].Overall, HW can be defined as a prolonged period of successive days (at least three) of abnormally high temperatures that exceed a specific thermal threshold of the local climate conditions of a specific region during the warm period (summer) of the year.The magnitude, duration, geographical extent, and severity of the HW are all taken into account [11,29,30].
There is no universally agreed upon definition of an HW, but there are many methods and indices used in the scientific literature to define HWs.HWs are sometimes defined using fixed temperature thresholds, which highlight the hot weather events that may have adverse health outcomes [18,30,31].However, fixed temperature thresholds are applicable to specified time periods at a particular site, and they may be applied to a region without considering the variations in latitude, climate, and season [32,33].Therefore, percentile thresholds are commonly used because their values vary with different climatic and geographical characteristics at each study site [34].For example, [33,35] recommended using the percentile threshold method to define HWs.Percentile-based threshold indices ranging from the 75th to the 99th are often used to determine temperature extremes [19,[36][37][38][39][40].Using a high percentile (99th) often does not cover all temperature extremes, whereas a low percentile (75th) may include a relatively large number of non-extreme events, so the middle percentile (90th) may be the more appropriate choice [41,42].
The increase in HWs and their aspects have global impacts on human mortality [43][44][45] and negative influences on the different social and economic sectors [33,46,47].The Middle East and North Africa (MENA) region is one of the most highly affected regions by global climate warming, which in turn leads to more extreme environmental conditions (hot and dry) than already exist [48,49].Most areas in MENA countries receive less than 300 mm of annual precipitation.This makes them extremely sensitive to changes in temperature and precipitation [49].Analysis of the long-term  temperature in the Eastern Mediterranean and Middle East region revealed that there have been increasing trends in the frequency and number of HWs since the 1970s [50].Moreover, the future projections of the duration and intensity of HWs in the MENA region during the 21st century, using two climate scenarios (RCP4.5 and RCP8.5), revealed that most of the region would experience unfavorable temperature conditions [51].Although the Kingdom of Saudi Arabia (KSA) already has a very hot and dry climate, some recent studies have indicated that there will be future significant upward warming trends, which will in turn contribute to an increase in the frequency of extreme weather events in KSA [52][53][54].To assess HWs in an integrated approach, [19] suggested a new HW structure comprising number (HMN), frequency (HWF), duration (HWD), and amplitude (HWA) elements, while [42] added the HW magnitude (HWM) aspect.
Several studies have been conducted to analyze the long-term trends of temperature and extreme warm events in KSA [54][55][56], mostly at stations and sometimes with spatial coverage.Nonetheless, there are no publications that focus on assessing the different aspects of HWs using percentile indices and interpreting their climatological and spatiotemporal variations.Therefore, this study aims to (i) assess the decadal changes in the different HW aspects over KSA during the last four decades from 1982 to 2021, (ii) interpret the spatiotemporal patterns, anomalies, and variations in the different HW aspects over the study area, (iii) investigate the time-latitudinal variations in the different HW aspects.

Study Area and Meteorological Data
The Kingdom of Saudi Arabia (KSA) territory (Figure 1) has desert climatic conditions (very hot and harsh) that vary greatly with geography and season [57][58][59].Some other studies have projected that the KSA will be exposed to increasing and continuous warming trends during the coming decades.These trends could lead to an increase in the frequency and intensity of extreme weather events [52,53].
Atmosphere 2023, 14, x FOR PEER REVIEW 3 of 21 warming trends during the coming decades.These trends could lead to an increase in the frequency and intensity of extreme weather events [52,53].The available meteorological stations in KSA are spatially and temporally dispersed, and they have a degree of uncertainty due to their proximity to urban communities, airports, and seaports.Additionally, desert areas such as the southeastern part of KSA, known as the Empty Quarter (647,500 km 2 ), do not have any meteorological stations [60,61].Therefore, relying on reanalysis datasets (such as ERA5) is highly beneficial for overcoming the absence, scarcity, and unavailability of both temporal and spatial meteorological observations for various weather and climate-related applications [11].Furthermore, [62] evaluated near-surface temperature (2 m) from ERA5 over several stations representing different climate zones in KSA.They found that ERA5 has a good agreement with the observed values of temperature, where the correlation coefficient ranged from 0.86 to 0.97 and the root mean square error ranged from 0.197 to 0.339 °C.Hence, the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth generation reanalysis climate (ERA5) dataset [63] is used to achieve the different purposes of this study.The required meteorological parameters to compute the five HW aspects are the 2 m temperatures (maximum (mx2t) and minimum (mn2t)) and total precipitation (tp).The ERA5 hourly data for mx2t, mn2t, and tp are freely available and downloaded on regular grids of 0.25° from the Copernicus Climate Change Service (C3S) Climate Date Store website provider (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels;accessed on 15 October 2022) over the KSA for the study period 1982-2021.ERA5 provides spatially and temporally homogeneous (long-term) climate variables that are generated from combinations of huge amounts of historical observations (in situ, satellite, etc.) into global models using an advanced data assimilation technique.The daily maximum for mx2t (dmxt2), daily minimum of mn2t (dmnt2), and the daily accumulation of tp (dtp) are computed from the obtained ERA5 hourly data using the climate data operator (CDO) tool.

Methods and Heatwave Characteristics
Various scientific research communities, experts, and scientists have developed several climate indices, including temperature indices, to serve different climate-related applications around the world.The Expert Team on Climate Change Detection and Indices (ETCCDI) developed 27 descriptive cores for climate change indices, including 16 cores for temperature indices (http://etccdi.pacificclimate.org/list_27_indices.shtml, accessed on The available meteorological stations in KSA are spatially and temporally dispersed, and they have a degree of uncertainty due to their proximity to urban communities, airports, and seaports.Additionally, desert areas such as the southeastern part of KSA, known as the Empty Quarter (647,500 km 2 ), do not have any meteorological stations [60,61].Therefore, relying on reanalysis datasets (such as ERA5) is highly beneficial for overcoming the absence, scarcity, and unavailability of both temporal and spatial meteorological observations for various weather and climate-related applications [11].Furthermore, Ref. [62] evaluated near-surface temperature (2 m) from ERA5 over several stations representing different climate zones in KSA.They found that ERA5 has a good agreement with the observed values of temperature, where the correlation coefficient ranged from 0.86 to 0.97 and the root mean square error ranged from 0.197 to 0.339 • C. Hence, the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth generation reanalysis climate (ERA5) dataset [63] is used to achieve the different purposes of this study.The required meteorological parameters to compute the five HW aspects are the 2 m temperatures (maximum (mx2t) and minimum (mn2t)) and total precipitation (tp).The ERA5 hourly data for mx2t, mn2t, and tp are freely available and downloaded on regular grids of 0.25 • from the Copernicus Climate Change Service (C3S) Climate Date Store website provider (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels;accessed on 15 October 2022) over the KSA for the study period 1982-2021.ERA5 provides spatially and temporally homogeneous (long-term) climate variables that are generated from combinations of huge amounts of historical observations (in situ, satellite, etc.) into global models using an advanced data assimilation technique.The daily maximum for mx2t (dmxt2), daily minimum of mn2t (dmnt2), and the daily accumulation of tp (dtp) are computed from the obtained ERA5 hourly data using the climate data operator (CDO) tool.

Methods and Heatwave Characteristics
Various scientific research communities, experts, and scientists have developed several climate indices, including temperature indices, to serve different climate-related applications around the world.The Expert Team on Climate Change Detection and Indices (ETCCDI) developed 27 descriptive cores for climate change indices, including 16 cores for temperature indices (http://etccdi.pacificclimate.org/list_27_indices.shtml, accessed on 16 July 2022), to monitor and detect the changes in climate extremes [64].The WMO committee for climatology, in consultation with representatives in different sectors (Water Resources and Hydrology, Health, and Agriculture and Food Security), established an Expert Team on Sector-Specific Climate Indices (ET-SCI) that identified and assessed an additional 34 sector-specific climate indices [42].The ET-SCI expert team built and developed an R-based program called ClimPACTv2 to calculate their defined sector-specific indices that require a homogeneous climate data series [65].The computed indices from the ClimPACTv2 software (ARC Centre of Excellence for Climate System Science, University of Melbourne, Melbourne, Australia) required the daily maximum air temperature (dmxt2), daily minimum air temperature (dmnt2), and daily total precipitation (dtp), which are obtained in this study from the ERA5 dataset.
The ClimPACTv2 software is implemented in this study to investigate the HW severity over KSA based on the calculated five HW aspects (HMN, HWF, HWD, HWM, and HWA).These HW aspects (Table 1) in ClimPACTv2 are computed using three HW indices, namely, the excess heat factor (EHF), the 90th percentile of daily maximum temperature (TX90), and the 90th percentile of daily minimum temperature (TN90), based on the reference climate period from 1982 to 2011.The ET-SCI expert team defined an HW in the ClimPACTv2 software as three or more successive days occurring in the summer (May-September in the Northern Hemisphere), in which either EHF > 0, TX > TX90, or TN > TN90.EHF (a measure of heatwaves and their intensity) was developed by [66] in the Australian Bureau of Meteorology, and recently augmented by [42,67].The sequential steps for calculating the EHF are described in detail in these aforementioned studies along with other recent studies [11,68].Total number of days that contribute to the annual HW events as identified by HWN.

HWD Heatwave Duration
Total number of days that contribute to the longest annual HW event identified by HWN.

HWM Heatwave Magnitude
Average of daily temperature across all annual HW events identified by HWN.

HWA Heatwave Amplitude
The highest daily temperature in the hottest annual HW event with highest HWM.
• C ( • C 2 for EHF) Moreover, the decadal sum of the identified five HW aspects for each HW index (EHF, TX90, and TN90) is aggregated from the yearly values within each decade (1982-1991, 1992-2001, 2002-2011, and 2012-2021), while the decadal anomalies for each HW aspect are computed as the decadal HW aspect subtracted from the average of its four decadal sums resulting from the same HW index.Also, the time-latitudinal cross-section (Hovmöller diagram) is used to study and interpret the annual behaviors and latitudinal variations in the five HW aspects during the study period (1982-2021) over KSA.The time-latitudinal cross-section is calculated as the sum of the HW aspect values at all longitudinal grid points within each latitudinal grid point for each year of the entire study period for each HW index.

Decadal Analysis of HWN
The spatial distribution of the decadal sum and anomaly of annual HWN based on TX90, TN90, and EHF indices throughout the last four decades (1982-1991, 1992-2001, 2002-2011, and 2012-2021) are shown in Figure 2. The first, third, and fifth rows in   Also, the range of the total occurred HWN (events) and their anomalies within each decade derived from TX90, TN90, and EHF indices are demonstrated in Table 2.It is detected that the decadal sum and anomaly of HWN increase spatially and temporally from the first to the fourth decade based on the three indices, as demonstrated by Table 2.In the first decade, most of KSA is invaded by less than 25 decadal HWN, with mostly negative anomaly, and the highest decadal sum and anomaly of HWN are found along the Red Sea coast, especially from TX90.The increases in decadal sum and anomaly of HWN are found on both spatial and temporal dimensions, starting from the second decade, with 28 events and the most negative anomaly (−15 to 5).The third decade has a decadal sum range of HWN from 12 to 42 events, with the most positive anomaly (−9 to 17), while it increases up to the highest decadal values in the fourth decade, with 17-60 events and considerable positive anomaly (5 to 34).Although the decadal sum and anomaly of HWN from the three indices show a marked gradual increase, their spatial distributions differ from one index to another.The highest spatial decadal sum and anomaly of HWN were produced from TX90 in the first two decades, from EHF in the third decade, and from TN90 in the last decade.Additionally, the southern part of KSA has the lowest decadal sum and anomaly of HWN, while their highest values are found along the Red Sea coast and in some different locations distributed in the north and middle of KSA.The results of increasing HWN coincide with the projected studies of HW events in North Africa, the Eastern Mediterranean, and the Middle East [11,52,53,69,70].
The time-latitudinal cross-section of annual HWN over KSA derived from TX90, TN90, and EHF indices is shown in Figure 3

Decadal Analysis of HWF
Figure 4 shows the decadal sum (first, third, and fifth rows) and decadal anomaly (second, fourth, and sixth rows) of annual HWF (all days that contribute to all HW events) during the selected four decades derived from the three (TX90, TN90, and EHF) indices.Also, Table 3 demonstrates the spatial extent range of the decadal sum and anomaly of HWF during each decade from the three indices.According to Figure 4 and Table 3, the lowest decadal sum of HWF (3:102 days) with all negative anomalies (−120: −17 days) was found in the first decade, followed by the second decade.The highest decadal sum of HWF (70:345 days) with all positive anomalies (22:225 days) was observed in the fourth decade, followed by the third decade.In terms of spatial distribution, TX90 has the highest decadal sum and anomaly of HWF, followed by EHF, in the first decade.EHF has the highest values during the other three decades, followed by TX90 in the second and third decades, and by TN90 in the last decade.Moreover, there is a considerable spatiotemporal gradual increase in the decadal sum and anomaly of HWF from the first to the fourth decade, as noticed from Figure 4 and Table 3.This remarkable temporal increase in the HWF agrees with several studies over the Middle East [52,53,71] and Egypt [11].
Additionally, the time-latitudinal cross-section of annual HWF has relatively the same behavior as the HWN, as shown in Figure 5.The period from 1982 to 1996 had the smallest annual HWF of less than 200 days, except for some latitudes where the three indices recorded high values (more than 600 days) in 1987, 1990, and 1991.From 1997 to 2008, the maximum annual HWF (more than 600 days) is detected in 1998 and 1999, followed by 2006 and 2007 (more than 300 days), while the lowest annual HWF (less than 300 days) is found in 2005.The highest annual HWF (more than 1000 days) across the study period is recorded during the last years from 2015 to 2021, followed by 2010 and 2012.Also, the time-latitudinal cross-section of annual HWF reveals that the annual HWF

Decadal Analysis of HWF
Figure 4 shows the decadal sum (first, third, and fifth rows) and decadal anomaly (second, fourth, and sixth rows) of annual HWF (all days that contribute to all HW events) during the selected four decades derived from the three (TX90, TN90, and EHF) indices.Also, Table 3 demonstrates the spatial extent range of the decadal sum and anomaly of HWF during each decade from the three indices.According to Figure 4 and Table 3, the lowest decadal sum of HWF (3:102 days) with all negative anomalies (−120: −17 days) was found in the first decade, followed by the second decade.The highest decadal sum of HWF (70:345 days) with all positive anomalies (22:225 days) was observed in the fourth decade, followed by the third decade.In terms of spatial distribution, TX90 has the highest decadal sum and anomaly of HWF, followed by EHF, in the first decade.EHF has the highest values during the other three decades, followed by TX90 in the second and third decades, and by TN90 in the last decade.Moreover, there is a considerable spatiotemporal gradual increase in the decadal sum and anomaly of HWF from the first to the fourth decade, as noticed from Figure 4 and Table 3.This remarkable temporal increase in the HWF agrees with several studies over the Middle East [52,53,71] and Egypt [11].
increases gradually from the beginning of the study period and reaches its maximum during the fourth decade after 2015.Additionally, the time-latitudinal cross-section of annual HWF has relatively the same behavior as the HWN, as shown in Figure 5.The period from 1982 to 1996 had the smallest annual HWF of less than 200 days, except for some latitudes where the three indices recorded high values (more than 600 days) in 1987, 1990, and 1991.From 1997 to 2008, the maximum annual HWF (more than 600 days) is detected in 1998 and 1999, followed by 2006 and 2007 (more than 300 days), while the lowest annual HWF (less than 300 days) is found in 2005.The highest annual HWF (more than 1000 days) across the study period is recorded during the last years from 2015 to 2021, followed by 2010 and 2012.Also, the time-latitudinal cross-section of annual HWF reveals that the annual HWF increases gradually from the beginning of the study period and reaches its maximum during the fourth decade after 2015.

Decadal Analysis of HWD
The spatial distribution of the decadal sum and anomaly of annual HWD during the selected four decades are shown in Figure 6.The upper two rows show the data from TX90, the middle two rows show the data from TN90, and the bottom two rows show the data from EHF.The ranges of the data are indicated in Table 4.The decadal sum and

Decadal Analysis of HWD
The spatial distribution of the decadal sum and anomaly of annual HWD during the selected four decades are shown in Figure 6.The upper two rows show the data from TX90, the middle two rows show the data from TN90, and the bottom two rows show the data from EHF.The ranges of the data are indicated in Table 4.The decadal sum and anomaly of HWD have almost the same behaviors as HWN and HWF.The lowest decadal sum of HWD (3:47 days) with negative decadal anomaly (−60:1 days) was found in the first decade.The HWD decadal sum gradually increased from the second decade to the highest value (34:145 days) with positive anomaly (4:85 days) in the fourth decade.This considerable increase in the HWD during the study period matches results identified by [53] over the Middle East and North Africa (MENA) and by [11] over Egypt.The largest decadal sum and anomaly of HWD over KSA area are obtained from TX90 followed by EHF in the first decade, from EHF followed by TX90 in the middle two decades, followed by TN90 in the fourth decade.In addition, the highest HWD is found along the Red Sea coast and the northern part of KSA, while the lowest HWD is found in the southern part of KSA.
Atmosphere 2023, 14, x FOR PEER REVIEW 10 of 21 [53] over the Middle East and North Africa (MENA) and by [11] over Egypt.The largest decadal sum and anomaly of HWD over KSA area are obtained from TX90 followed by EHF in the first decade, from EHF followed by TX90 in the middle two decades, followed by TN90 in the fourth decade.In addition, the highest HWD is found along the Red Sea coast and the northern part of KSA, while the lowest HWD is found in the southern part of KSA.The time-latitudinal cross-section of the annual HWD calculated from the three indices (TX90, TN90, and EHF) is shown in Figure 7, which has relatively the same HWN and HWF behaviors.The annual HWD during the period from 1982 to 1996 (often less than 200 days) is lower than in the rest of the study period, and the maximum values (more than 400 days) are detected after 2015.This implies that there is a considerable temporal gradual increase in annual HWD from the first to the fourth decade, where the highest values across the study period are recorded in 1987, 1991, 1998, 2003, 2010, 2012, and   The time-latitudinal cross-section of the annual HWD calculated from the three indices (TX90, TN90, and EHF) is shown in Figure 7, which has relatively the same HWN and HWF behaviors.The annual HWD during the period from 1982 to 1996 (often less than 200 days) is lower than in the rest of the study period, and the maximum values (more than 400 days) are detected after 2015.This implies that there is a considerable temporal gradual increase in annual HWD from the first to the fourth decade, where the highest values across the study period are recorded in 1987,1991,1998,2003,2010,2012, and from 2015 to 2021.

Decadal Analysis of HWM
The average daily temperatures (from TX90 or TN90) in • C, or the average daily EHF values in • C 2 (from EHF), over all days contributed to all annual HW events express the annual HWM. Figure 8 shows the spatial patterns of the decadal sum (first, third, and fifth rows) and anomaly (second, fourth, and sixth rows) of HWM during the chosen four decades over KSA.The ranges of the HWM decadal sum and anomaly obtained from the three indices (TX90, TN90, and EHF) over KSA for each decade are illustrated in Table 5.

Decadal Analysis of HWM
The average daily temperatures (from TX90 or TN90) in °C, or the average daily EHF values in °C2 (from EHF), over all days contributed to all annual HW events express the annual HWM. Figure 8 shows the spatial patterns of the decadal sum (first, third, and fifth rows) and anomaly (second, fourth, and sixth rows) of HWM during the chosen four decades over KSA.The ranges of the HWM decadal sum and anomaly obtained from the three indices (TX90, TN90, and EHF) over KSA for each decade are illustrated in Table 5.It is noticed that the decadal sum and anomaly of HWM derived from TX90 is greater than that derived from TN90 and EHF (producing the lowest values) during all decades, which agrees with the results obtained by [11] over Egypt.The lowest decadal sum of HWM with predominantly negative decadal anomalies from the three indices is found during the first decade followed by the second decade, especially over the southern and southwestern parts of KSA.The decadal sum of HWM from TX90 (first row in Figure 8) over most of KSA during the first two decades is more than 200 • C, with mostly negative (less than −100 • C) anomalies (second row in Figure 8), while it is generally greater than 350 • C, with mostly positive (more than 50 • C) anomalies, in the last two decades.For the TN90 index, the decadal sum of HWM (third row in Figure 8) ranges from 25 to 250 • C, with mostly negative anomalies (less than −20 • C), in the first two decades.In the last two decades, the decadal sum ranges from 200 to 300 • C, with mostly positive anomalies (more than 30 • C).For the EHF index, the decadal sum of HWM (fifth row in Figure 8) is less than 50 • C 2 over most parts of KSA, especially in the first two decades.It is less than 10 • C 2 over the southern parts in the last two decades, while it ranges from 20 to 70 • C 2 over the northern parts in the last two decades.The decadal anomaly of HWM derived from EHF is mostly negative (−20:8 • C 2 ) in the first two decades and mostly positive (−6:30 • C 2 ) in the last two decades, where it varies greatly in the northern and northwest parts of KSA.Moreover, the decadal HWM increases gradually over KSA from the first to the last decade based on the three utilized indices.
Figure 9 shows the annual HWM time-latitudinal cross-section, where the highest, moderate, and lowest values are derived from TX90 (exceeds 1800 • C), TN90 (less than 1300 • C), and EHF (less than 250 • C 2 ), respectively.The time-latitudinal distribution of annual HWM behaves in a similar way for HWN, HWF, and HWD aspects, where the lowest values appear during the period 1982-1996.After 1996, the annual HWM increases gradually over time, and the maximum values are found during the period from 2015 to 2021, especially over the central and northern parts of KSA.Moreover, the maximum values of time-latitudinal HWM before 1996 were found in 1987 and 1988.The maximum values appeared in 1998, 2003, 2010, and 2012 during the period 1997-2014.However, the highest maximum values were detected across the period 2015-2021 based on the three indices.These findings are consistent with the results of [11,53], as they show that HWM is one of the crucial factors exacerbating HWs in the Middle East and North Africa.

Decadal Analysis of HWA
The highest daily temperature in °C (from TX90 or TN90), or highest daily EHF value in °C2 (from EHF), appeared across all annual HW events, indicating the HWA.The spatial distribution of the decadal sum and anomaly of the annual HWA over KSA is shown in Figure 10, as derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows) during the four study decades (columns from left to right).Also, Table 6 shows the available ranges of the HWA decadal sum and anomaly over KSA during each decade based on TX90, TN90, and EHF indices.It is found that the spatial patterns of decadal sum and anomaly of the annual HWA are remarkably similar to the corresponding spatial patterns of the HWM obtained from the three indices.The smallest decadal sum of HWA with the maximum negative anomaly was found in the first decade, followed by the second decade.This was especially true in the northern part of KSA from the TX90 index (sum < 200 °C; anomaly < −60 °C) and the EHF index (sum < 20 °C2 ; anomaly < 0 °C2 ).It was also found in several zones in the southwestern, northeastern, and central parts of KSA from the TN90 index (sum < 100 °C; anomaly < −60 °C).The largest HWA decadal sums with maximum positive anomaly over most of KSA from TX90 (sum > 300 °C; anomaly > 40 °C), TN90 (sum > 250 °C; anomaly > 50 °C), and EHF (sum > 40 °C2 ; anomaly > 15 °C2 ) are found in the fourth decade, followed by the third decade.These obtained results are in agreement with the findings of projected changes in HWA studied by [11,69,70].In addition, the largest decadal sum and anomaly of HWA during the four decades over KSA are obtained from TX90, followed by TN90 and EHF, respectively.

Decadal Analysis of HWA
The highest daily temperature in • C (from TX90 or TN90), or highest daily EHF value in • C 2 (from EHF), appeared across all annual HW events, indicating the HWA.The spatial distribution of the decadal sum and anomaly of the annual HWA over KSA is shown in Figure 10, as derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows) during the four study decades (columns from left to right).Also, Table 6 shows the available ranges of the HWA decadal sum and anomaly over KSA during each decade based on TX90, TN90, and EHF indices.It is found that the spatial patterns of decadal sum and anomaly of the annual HWA are remarkably similar to the corresponding spatial patterns of the HWM obtained from the three indices.The smallest decadal sum of HWA with the maximum negative anomaly was found in the first decade, followed by the second decade.This was especially true in the northern part of KSA from the TX90 index (sum < 200 • C; anomaly < −60 • C) and the EHF index (sum < 20 • C 2 ; anomaly < 0 • C 2 ).It was also found in several zones in the southwestern, northeastern, and central parts of KSA from the TN90 index (sum < 100 • C; anomaly < −60 • C).The largest HWA decadal sums with maximum positive anomaly over most of KSA from TX90 (sum > 300 • C; anomaly > 40 • C), TN90 (sum > 250 • C; anomaly > 50 • C), and EHF (sum > 40 • C 2 ; anomaly > 15 • C 2 ) are found in the fourth decade, followed by the third decade.These obtained results are in agreement with the findings of projected changes in HWA studied by [11,69,70].In addition, the largest decadal sum and anomaly of HWA during the four decades over KSA are obtained from TX90, followed by TN90 and EHF, respectively.The time-latitudinal cross-section of the annual HWA over KSA during the study period (1982-2021) is shown in Figure 11, where it mostly has the same behavior of HWM.The time-latitudinal cross-section of the annual HWA over KSA during the study period (1982-2021) is shown in Figure 11, where it mostly has the same behavior of HWM.The highest time-latitudinal HWA values are obtained from TX90, while the moderate and the lowest values are obtained from TN90 and EHF, respectively.Also, the lowest time-   Generally, the observed increase in the different HW aspects coincides with the behavior and trend of temperature over KSA and HW over the surrounding region, which is confirmed by several previous studies, e.g., [54,[72][73][74][75][76][77].

Limitations of this Study
The 90th percentile indices employed in this study may perform well in determining heatwaves across KSA, but not all estimated heatwaves may have an impact on different sectors, such as health, agriculture, etc.
The evaluation of ERA5 reanalysis dataset accuracy was restricted to the available stations only due to the unavailability of measured data at more ground stations in KSA; however, the ERA5 dataset showed good accuracy at those few available stations.In addition, the availability of measured data may greatly contribute to a more accurate interpretation and analysis of HW aspects in KSA.

Conclusions
In this study, we assessed the decadal changes, spatio-temporal distribution, anomalies, and variations in heatwave (HW) aspects over Saudi Arabia (KSA) during the last four decades .These HW aspects were computed by the ClimPACT2 (R-based) software based on three threshold indices (TX90, TN90, and EHF), using the prepared daily gridded (0.25 • × 0.25 • ) maximum (TX) and minimum (TN) temperatures from the ERA5 reanalysis dataset.These HW aspects are number (HWN), frequency (HWF), duration (HWD), magnitude (HWM), and amplitude (HWA).The main findings of this study can be summarized as follows: -There is a remarkable congruence in the behaviors and patterns of the five HW aspects, whether they were derived from the TX90, TN90, or EHF index, although they produced relatively different estimates for these aspects.The highest maximum values of all HW aspects are found during the period 2015-2021, while the lowest values are detected in the period 1984-1986, followed by 1992 and 1993.
On the other hand, future work may focus on assessing the impact of HWs on various KSA sectors, as well as studying thermal comfort and heat stress indices.

Figure 1 .
Figure 1.Topographical map and features of the KSA territory inside the solid bold polygon.

Figure 1 .
Figure 1.Topographical map and features of the KSA territory inside the solid bold polygon.

Figure 2
Figure 2 illustrate the decadal sum of annual HWN, while the second, fourth, and sixth rows show the decadal anomaly of annual HWN derived from TX90, TN90, and EHF indices, respectively.

Figure 2 .
Figure 2. Decadal sum and anomaly of heatwave number derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows).

Figure 2 .
Figure 2. Decadal sum and anomaly of heatwave number derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows).
. The years from 1982 to 1996 have the lowest annual HWN, with maximum time-latitudinal annual HWN of less than 140 events, especially in 1987, 1990, and 1991 from the three indices.The maximum time-latitudinal annual HWN (more than 200 events) during the period from 1997 to 2008 is detected in 1998, followed by 1999, 2000, 2006, and 2007 (more than 100 events).The last period from 2009 to 2021 has the highest time-latitudinal annual HWN (more than 200 events), particularly in 2010, 2015-2019, and 2021, with maximum values (more than 250 events) in 2017, 2019, and 2021.

Figure 4 .
Figure 4. Decadal sum and anomaly of heatwave frequency derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows).

Figure 4 .
Figure 4. Decadal sum and anomaly of heatwave frequency derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows).

Figure 6 .
Figure 6.Decadal sum and anomaly of heatwave duration derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows).

Figure 6 .
Figure 6.Decadal sum and anomaly of heatwave duration derived from TX90 (first and second rows), TN90 (third and fourth rows), and EHF (fifth and sixth rows).
latitudinal annual HWA values occur in the first decade, with annual increasing trend behavior up to the highest values in the last (fourth) decade, based on the three indices.According to the three indices, the largest time-latitudinal HWM values were detected in 1987 and 1988 during the first decade, in 1998 during the second decade, and in 2003, 2010, and 2012 during the third decade.The highest values appeared in the period 2015-2021 during the last decade.

-
The decadal sum and anomaly of the five HW aspects over the entire KSA territory increase gradually from the first decade(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991) to the highest values in the last decade (2012-2021), based on the three indices.-The highest values of these HW aspects are found mostly in the northern and western parts of KSA, particularly along the Red Sea coast, while the lowest values are found in the southern part of KSA.-The maximum decadal sum with mostly positive anomaly of HWN (42 events), HWF (255 days), and HWD (145 days) occurred in the last decade, as derived from the three indices.-The maximum decadal sum of HWM (HWA) ranges from 175 to 463 (from 189 to 471) • C based on TX90 and TN90, while it ranges from 7 to 58 (from 15 to 185) • C 2 based on EHF.-The time-latitudinal cross-section analysis revealed that the maximum values of the five HW aspects are found in 1983, 1987, and 1991 (first decade), 1998 (second decade), and 2003, 2006, 2007, and 2010 (third decade), as well as 2012 and 2015-2021 (last decade).Also, the influence of HW aspects from all indices is most pronounced in the northern part of the KSA in most decades, especially in the last decade.- and editing, M.M., A.A. and H.A.B.All authors have read and agreed to the published version of the manuscript.

Table 2 .
The decadal sum and anomaly ranges for heatwave number from TX90, TN90, and EHF indices.

Table 3 .
The decadal sum and anomaly ranges for heatwave frequency from TX90, TN90, and EHF indices.

Table 3 .
The decadal sum and anomaly ranges for heatwave frequency from TX90, TN90, and EHF indices.

Table 4 .
The decadal sum and anomaly ranges for heatwave duration from TX90, TN90, and EHF indices.

Table 4 .
The decadal sum and anomaly ranges for heatwave duration from TX90, TN90, and EHF indices.
The highest time-latitudinal HWA values are obtained from TX90, while the moderate and the lowest values are obtained from TN90 and EHF, respectively.Also, the lowest time-latitudinal annual HWA values occur in the first decade, with annual increasing trend behavior up to the highest values in the last (fourth) decade, based on the three indices.According to the three indices, the largest time-latitudinal HWM values were detected in 1987 and 1988 during the first decade, in 1998 during the second decade, and in 2003, 2010, and 2012 during the third decade.The highest values appeared in the period 2015-2021 during the last decade.