Initial Spatio-Temporal Assessment of Aridity Dynamics in North Macedonia (1991–2020)

1. Introduction

In recent decades, climate variability and climate change have been central topics in numerous scientific studies. Global climate change, driven by both natural processes and anthropogenic influences, represents a major environmental challenge of the twenty-first century [1,2,3,4,5]. Increased attention from governments, the scientific community, and the general public underscores the urgency of understanding these changes and their effects [6,7]. Temperature fluctuations, the rate of change, and associated phenomena such as rising sea levels, floods, droughts, and aridity underscore the broader implications of climate change [8]. Seasonal analyses reveal a consistent warming trend in near-surface air temperatures, with the most pronounced increase occurring during winter [5].

Aridity plays a crucial role in understanding climatic conditions and natural processes, serving as both an indicator of plant growth and a determinant of landscape characteristics. Evaluation of the aridity can be managed through the widely used specific climatic parameters. Parameters in which there is a defined relationship between air temperature and precipitation are known as “aridity indices” and represent a useful tool in their quantification [9]. However, formulation of aridity indices is not straightforward due to the nature of aridity. Use of mean annual precipitation as an index of aridity ignores the importance of air temperature and evaporation. Aridity indexes that have gained widespread acceptance directly or indirectly take into account all three factors [10].

To quantify aridity, numerous climatic indices have been formulated, generally calculated as a ratio between some specific cli¬matic indicators. For example, the Lang Rain Factor [11], De Martonne Aridity Index [12], Emberger Aridity Index [13], Gračanin’s rain factor for the vegetation period (1950) [14], Seljanin’s hydrothermal coefficient (1958) [15], Erinç aridity index (Erinç, 1965) [16] and Pinna Combinative Aridity Index (Zambakas, 1992) [17] are calculated as the ratio of precipitation and air temperature. Since these indices are derived from different climatic variables (calculated at monthly, annual scale, vegetation period), they provide heterogenous perspectives in aridity analysis. Also, there are some different approaches in determining aridity such as: Thornthwaite aridity index (1948) [18], Budyko’s equation (1951) [19], UNESCO aridity index (UNESCO, 1979) [20] and FAO aridity index (FAO, 1989) [21]. Thornthwaite (1948) [18] introduced an aridity index as the ratio between annual precipitation and potential evapotranspiration. A few years later, Budyko (1951) [19] offered a new approach by considering the heat and water balance equations of the Earth’s surface. The Radiational Index of Dryness (RID) was the ratio of the mean annual net radiation (i.e. the radiation balance) to the product of the mean annual precipitation times the latent heat of vaporization for water. The UNESCO aridity index is based on the ratio of average annual precipitation to potential evapotranspiration. This index, originally developed by UNESCO (1979) [20] to map drylands, remains a widely used tool for assessing and projecting global aridity and its impacts. The FAO Aridity Index (1989) [21] measures a region's dryness by calculating the ratio of annual precipitation to potential evapotranspiration. This index is a critical tool for assessing moisture availability, managing land resources for agriculture, and understanding the long-term impacts of climate change on water resources and drylands.

Furthermore, some other indices, such as the Johansson Continentality Index and the Kerner Oceanity Index, differentiate continental and maritime climate [22,23]. Recent advances in the archiving and collection of data have allowed worldwide monitoring of aridity as opposed to hindsight assessments. The most commonly used satellite measure is the Normalized Difference Vegetation Index (NDVI). The assessment of aridity changes using remote sensing like NDVI measures vegetation greenness, which declines in dry conditions. However, this requires images of medium-high spatial resolution over a long period of time and high-quality data derived from Landsat [10,24].

Since the 1950s, aridity has increased across most of the Earth’s surface, a trend intensified by the ongoing effects of global warming induced by anthropogenic activities [25]. Sardans et al. (2024) [26] concluded that from 1960 to 2023, 27.9% of the global land surface became significantly more arid, while 20.5% became significantly less arid. This indicates a shift towards drier climates, with humid, semi-humid, and semi-arid areas decreasing by 8.51, 1.45, and 0.53 million-km², respectively, and arid and hyper-arid areas increasing by 6.34 and 4.18 million-km², respectively. Especially vulnerable regions on global scale according to this study are Europe, North Africa and the Middle East. The spatio-temporal dynamics and variability of the aridity was investigated in details on regional scale: Greece [27,28,29], Italy [30,31], the Iberian Peninsula [32,33], Turkey [23.34], Algeria [35], the Middle East [36], Lebanon [37], Iran [38,39] and Iraq [40]. Thornthwaite (1948) [18] introduced an aridity index as the ratio between annual precipitation and potential evapotranspiration. Other indices, such as the Johansson continentality index and the Kerner oceanity index, differentiate between continental and oceanic climates [22,23].

Additionally, Mărgărit-Mircea (2021) [41] used the Johansson Index of Continentality, the Kerner Index of Oceanity, the De Martonne Aridity Index, the Pinna Combinative Index and the UNEP Aridity Index. This study determined that the Southeastern Europe, the Balkan Peninsula, the Iberian Peninsula, the Apennine Peninsula, and Mediterranean islands consistently exhibit a dry climate according to all indices. Also, Beštáková et al. (2023) [42] assessed the spatial distribution of arid and humid regions in Europe using the Aridity Index (1950–2019) identifying the Balkan Peninsula as predominantly arid. Numerous authors dealt with aridity occurrence in the Balkan countries, for instance in Serbia (e.g., [43,44,45,46]), Montenegro [47], Croatia [48,49], Romania (e.g. [50,51,52]), Bulgaria [53].

However, it is difficult to identify similar and robust aridity research on national scale. Aridity in North Macedonia has been examined by using FAO Aridity Index in Agro-Ecological Atlas of North Macedonia [54]. In this study, authors revealed the classes of semi-arid and dry sub-humid climate. These two classes delineate most of the agriculture areas in the central and southeast part of the country, which covers the driest areas. On the other hand, the humid class covers also some important agriculture areas in the western part of the country with more humid climatic condition. Recent studies indicate increasing of the future aridity conditions (RCP8.5 scenario) in the Southeastern Europe, including territory of North Macedonia. The Balkan Peninsula, including North Macedonia, has experienced a moderate increase in extreme climate events, particularly since the 1980s [55]. Additionally, aridity trends are interpreted within the broader context of paleoclimatic phenomena, particularly those occurring during the Middle and Late Pleistocene. Carvalho et al. (2022) [56] suggest strong reduction of Aridity Index for South and Southeastern Europe (2081-2100). Furthermore, [57,58] observed semy-dry conditions (2071-2100) according to De Martonne classes in central, south and southeastern parts of North Macedonia. The lack of adequate recent studies and above-mentioned facts justify the importance of reseach the spatio-temporal aridity patterns in North Macedonia.

For the purpose of this study, the Johansson Continentality Index, the De Martonne Aridity Index and the Pinna Combinative Index were selected as the most suitable tools for analyzing spatio-temporal conditions in North Macedonia. These climate indices were calculated using air temperature and precipitation data from 13 meteorological stations over a 30-year standard climatological interval (1991-2020) providing insight into the aridity dynamics on national scale. By examining climate data spanning several decades, this research aims to provide valuable insights into the region’s recent climate variability. Additionally, the objectives of this study are: (1) to analyze and discuss the trend characteristics of selected climate indices in details using Mann-Kendall non-parametric test; (2) to quantify the significance of changes on monthly, seasonal and annual scale; (3) to provide insight into spatial changes of indices using GIS tools and suitable interpolation technique, (4) to correlate our findings with the North Atlantic Oscillation Index (NAOI) and Southern Oscillation Index (SOI). By evaluating long-term trends and identifying regional climate shifts, this research contributes to the broader understanding of climate variability and its potential impacts on natural ecosystems and human societies.

2. Materials and Methods

2.1. Study area and data

North Macedonia, covering approximately 25,713 km² in the southern Balkan Peninsula, is characterized by a predominantly mountainous landscape, accounting for nearly 79% of its territory. In contrast, lowland plains constitute approximately 19%, while lakes and reservoirs occupy the remaining 2% of the total land area (Figure 1). The region exhibits pronounced geomorphological variability, with frequent alternations between mountainous terrain and valley systems, resulting in a relatively high average slope of 15.4° [59]. The mountain ranges of North Macedonia are arranged in sequences, forming effective barriers that reduce the inflow of moist air masses, particularly from the west. This orographic influence increases aridity in the central part of the country. The country's climatic conditions are shaped by the interplay of Mediterranean, continental, and orographic influences, contributing to significant spatial and temporal variations in air temperature and precipitation patterns [60]. The semiarid nature of the region is reflected in an average annual precipitation of approximately 500–700 mm. For the purposes of this study, observational data from 13 meteorological stations operated by the Hydrometeorological Service of North Macedonia (UHMR) [61] were used (see Table 1). However, ongoing climatic shifts have intensified extreme weather events, with an observable increase in the frequency and magnitude of convective storms and heavy precipitation episodes over recent decades.

Based on the updated Köppen–Geiger classification [62], North Macedonia exhibits a complex climatic pattern resulting from the interaction of topography, continental influences, and Mediterranean air masses (see Table 2). Lowland regions and river valleys, including the Vardar and Strumica basins, are characterized by a Cfa climate, corresponding to a humid subtropical regime with hot, dry summers and mild, wet winters, and mean annual air temperatures around 11–12 °C. Higher western and northwestern areas, such as those surrounding Bitola, fall under the Cfb subtype, representing a temperate oceanic climate with moderate summers, cool winters, and relatively uniform precipitation throughout the year, with mean annual air temperatures near 10–11 °C. Northern and eastern parts of the country are dominated by the Dfb subtype, a humid continental climate with cold winters, warm summers, and pronounced seasonal air temperature contrasts, with mean annual air temperatures around 8 °C. This distribution reflects both altitudinal and latitudinal gradients, underscoring the spatial variability of climate across North Macedonia, and aligns with observed trends of regional warming and precipitation variability over recent decades.

2.2. Pettitt’s homogeneity test

To evaluate the homogeneity of the climatic time series and identify potential abrupt shifts, we employed Pettitt's nonparametric change-point detection method [63]. This approach is particularly valuable for pinpointing the onset of statistically meaningful alterations in long-term records, such as those of air temperature and precipitation, where gradual trends may be confounded by sudden discontinuities.

Pettitt's test operates under the premise that a change point divides the series of n observations into two segments: the initial portion (up to time t) governed by cumulative distribution F1, and the subsequent segment by F2. A true shift implies that F1 ≠ F2, manifesting as a non-random deviation in the central tendency at some point τ. The null hypothesis (H₀) posits a single, uniform distribution across the entire series, indicative of homogeneity. The alternative (H1) posits a structural break at τ, driven by an exogenous factor altering the mean or median.

The test statistic for a given split t is the Mann-Whitney-type U_t,n, defined as:

$$U_{t,n}={\sum }_{i=1}^t{\sum }_{j=i+1}^{n}sgn\left(xi-xj\right)$$

where sgn denotes the sign function (+1 for positive differences, -1 for negative, and 0 otherwise). The overall test statistic K is then the maximum absolute deviation:

$$K_t=\underset{1\le t<n\left|U_{\begin{array}{c}t,n\end{array}}\right|}{\max }$$

under H₀, K follows an asymptotic distribution approximated by:

$$\rho =e^{-K/(n^2+n^3)}$$

where ρ represents the two-tailed significance probability. Rejection of H₀ occurs when p falls below a chosen α level (e.g., 0.05), with the corresponding t identifying the change point. All datasets underwent rigorous quality assurance by the UHMR [61]. We attribute any identified discontinuities to climatic signals rather than observational artifacts.

2.3. Johansson Continentality Index (JCI)

The Johansson Continentality Index (JCI) serves as a quantitative measure of the degree to which continental influences shape regional climatic conditions. It has been widely adopted in climatological and biogeographical research to characterize the transition between maritime and continental climates [22,64]. The index is determined using the following expression:

$$JCI=1.7\times \left({\displaystyle \frac{A}{\sin \phi }}\right)-20.4$$

where: A represents the annual air temperature amplitude, calculated as the difference between the mean air temperature of the warmest and coldest months (°C); and φ denotes the geographical latitude of the location, expressed in radians. This formulation integrates both thermal variability and latitudinal position, enabling the classification of climatic regimes along a continuum from oceanic to strongly continental types. The corresponding intervals of JCI values and their interpretative classes are summarized in Table 3.

2.4. The De Martonne aridity index

The De Martonne Aridity Index (I_DM) is a well-established and extensively utilized metric in applied climatology for assessing aridity and humidity conditions [66,67]. Despite being one of the earliest aridity indices, it remains a fundamental tool in climatic classification and is widely applied in contemporary research to delineate dry and humid regions across diverse geographical settings [22,35,68]. The annual I_DM is defined as follows:

$$I_{DM}=P÷\left(T+10\right)$$

where, P represents the total annual precipitation (mm), and T denotes the mean annual air temperature (°C). The inclusion of a constant term in the denominator ensures numerical stability and prevents division by low-temperature values. The index can also be computed at shorter temporal scales, such as monthly and seasonal intervals. The monthly I_DM is given by:

$${Im}_{\left(DM\right)}=\left(12\times P_m\right)÷\left(T_m+10\right)$$

where P_m corresponds to the monthly precipitation (mm), and Tm represents the mean monthly air temperature (°C). For seasonal assessments, the index is expressed as:

$${Is}_{\left(DM\right)}=\left(4\times P_s\right)÷\left(T_s+10\right)$$

where P_s and T_s refer to the seasonal precipitation (mm) and mean seasonal air temperature (°C), respectively. De Martonne's classification system categorizes climatic conditions based on I_DM values, providing a systematic framework for distinguishing between arid, semi-arid, humid, and extremely humid regions. The classification scheme is presented in Table 4.

De Martonne developed a climate classification based on the values of the I_DM, as shown in Table 4. This classification illustrates that an increase in I_DM values, assuming a constant air temperature, corresponds to an increase in precipitation, and vice versa. Essentially, higher I_DM values reflect more humid conditions, while lower values indicate drier climates. This relationship emphasizes the role of precipitation in shaping the moisture balance of a region, with I_DM acting as a key indicator of climate classification.

2.5. The Pinna combinative index

Pinna’s combinative or multiannual aridity index [7,17,22,23,69] represents a climatological metric designed to integrate both long-term and seasonal hydrothermal characteristics. It is formally expressed as:

$$I_p=\left({\displaystyle \frac{1}{2}}\right)*\left({\displaystyle \frac{P_{ma}}{T_{ma}+10}}+{\displaystyle \frac{12*P_d^{\text{'}}}{T_d^{\text{'}}+10}}\right)$$

where:

P_ma denotes the total annual precipitation (mm),

T_ma represents the mean annual air temperature (°C),

P’_d is the total precipitation during the driest month (mm),

T’_d refers to the mean air temperature of the driest month (°C).

According to the interpretative framework proposed by Baltas (2010) [70], values of I_p < 10 signify arid climatic regimes. When I_p ranges between 10 and 20, the climate is categorized as semiarid Mediterranean, typically associated with sclerophyllous vegetation characteristic of the formal Mediterranean biome. In contrast, values of I_p > 20 denote humid conditions (Table 5).

Unlike the I_DM, which allows for intra-annual temporal resolution, Pinna’s index is specifically calibrated for multiannual analyses. As such, it serves as a robust indicator for identifying climatological dryness during the most precipitation-deficient period of the year—an aspect of critical relevance in the context of agricultural planning and irrigation management.

2.6. Mann – Kendall trend test

Another widely used method for detecting temporal trends in hydrometeorological time series is the Mann–Kendall (MK) test, a non-parametric statistical approach commonly applied to air temperature, precipitation, and discharge records. In this study, it was employed to evaluate both the direction and the significance of monotonic trends, while the magnitude of the detected trends was estimated using Sen’s slope estimator [71]. The test assesses the relative ordering of data values by examining the number of positive and negative differences between subsequent observations. If a trend is present, the cumulative sign values exhibit a consistent tendency to either increase or decrease. The null hypothesis (H₀) posits the absence of a trend, whereas the alternative hypothesis (H₁) assumes the existence of a monotonic trend.

$$S={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^{n}sgn\left(x_i-x_j\right)$$

where the sign function sgn(x_j−x_i) is defined as

$$sgn\left(x\right)=\left\{\begin{array}{c}1ifx>x_i\\ 0ifx=x_i\\ -1ifx<x_i\end{array}\right.$$

The variance of S, denoted as VAR (S), is given by:

$$VAR\left(S\right)={\displaystyle \frac{1}{18}}\left[n\left(n-1\right)\left(2n+5\right)-{\sum }_{p-1}^g{t}_p\left(t_p-1\right)\left(2t_p+5\right)\right]$$

where n is the sample size, g is the number of tied groups, and t_p represents the number of observations in the ptℎ group. The standardized test statistic Z_MK is calculated as:

$$Z_{MK}=\left\{\begin{array}{c}{\displaystyle \frac{S-1}{\sqrt{VAR\left(S\right)}}}ifS>0\\ 0ifS=0\\ {\displaystyle \frac{S+1}{\sqrt{VAR\left(S\right)}}}ifS<0\end{array}\right.$$

A positive (negative) value of Z_MK indicates an increasing (decreasing) trend over time. For a given significance level α₁ where 0<α<0.5, and Z₁−α is the (1−α)-quantile of the standard normal distribution, the null hypothesis H₀ is rejected in favor of H1 when:

$$\begin{gathered} Z_{MK}=\ge Z_{1-\alpha } \\ {Z}_{MK}=\le Z_{1-\alpha } \end{gathered}$$

2.7. GIS software

To assess aridity in North Macedonia using the JCI, I_DM, and I_P indices, a custom C# software tool was developed to compute both indices annually for 13 meteorological stations, classifying their climate categories [72]. For spatial analysis, QGIS (v.3.38), an open-source geographic information system, was employed to map and interpolate aridity indices using the ordinary kriging method via its integration with SAGA GIS (v. 9.9.1) tools. Ordinary kriging was chosen for its ability to incorporate autocorrelation, considering both spatial arrangements and distances between points to minimize estimation errors, generating surface rasters with a spherical semivariogram model [73,74]. This approach, combining C# for index calculations with QGIS (v. 3.38) and SAGA (v. 9.9) for spatial interpolation and visualization, provides a robust framework for analyzing aridity trends across North Macedonia (Figure 2).

Pearson’s correlation coefficients [75] and associated p-values were computed using Python's scipy.stats.pearsonr function, which calculates r_p as:

$$r_p={\displaystyle \frac{{\sum }_{i=1}^n(x_i-\overline{x})(y_i-\overline{y)}}{\sqrt{{\sum }_{i=1}^n}(x_i-\overline{x}{)}^2\sqrt{{\sum }_{i=1}^n}(y_i-\overline{y}{)}^2}}$$

where xi represents I_DM values and y_i represents NAOI or SOI values, with means x̄ and ȳ. Statistical significance was tested at the 5% level (p < 0.05). All calculations were carried out in Python 3 using the SciPy package, specifically the Pearson’s function from the scipy.stats module.

3. Results and discussion

3.1. Pettitt’s homogeneity test

Before all calculation procedures the Pettitt’s homogeneity test was performed (Figure 3). Application of Pettitt’s nonparametric test to annual climate series (1991–2020) from 13 meteorological stations revealed contrasting change-point behavior for air temperature and precipitation. Annual mean air temperature (T_mean) exhibited a highly coherent response, with a statistically significant change point detected in 2006 at all stations. Test statistics (K_t) ranged from 160 to 204, with p-values between 0.00013 and 0.00406, indicating rejection of the homogeneity hypothesis at the 0.05 level. Positive U_t,n discontinuities indicate an abrupt post-2006 warming, with estimated increases of 0.8–1.5 °C relative to the pre-change period. The uniform timing and magnitude of these shifts suggest a dominant regional-scale forcing, consistent with the accelerated warming observed in the Balkan region since the early 2000s.

In contrast, annual total precipitation (P_total) showed pronounced spatial heterogeneity. Significant change points were identified at four stations (Mavrovo, Demir Kapija, Berovo, and marginally Prilep), all occurring in 2001, with Kt values between 117 and 149 and p-values below 0.05, except for the borderline case at Prilep. Negative U_t,n discontinuities indicate post-change decreases of approximately 150–300 mm yr⁻¹. No statistically significant shifts were detected at the remaining nine stations (p = 0.067–0.522). These results highlight a clear divergence between temperature and precipitation responses, with coherent regional warming contrasted by spatially variable and locally controlled precipitation changes.

3.2. Air temperature, precipitation, and aridity classification

To evaluate the spatial variation of aridity across North Macedonia, mean annual air temperature (MAT) and mean annual precipitation (MAP) were analyzed for the 1991–2020 reference period, providing critical insights into regional climatic patterns (Figure 4). These parameters underpin the calculation of the I_DM. The I_DM quantifies aridity, with higher values indicating more humid conditions and lower values denoting aridity. Data from 13 meteorological stations, spanning diverse topographic settings from low-elevation valleys to high-altitude mountainous regions, were used to assess spatial variability in these parameters.

The MAP values exhibit pronounced spatial heterogeneity, driven by topographic and geographic factors. The lowest MAP, 472.1 mm, is recorded at Štip, a station in the eastern lowlands at approximately 300 m elevation. This low precipitation reflects the region’s continental climate and rain-shadow effect, as it lies leeward of the western mountains. In contrast, the highest MAP, 1,058.5 mm, is observed at Lazaropole, located in the western mountainous region at <1,330 m elevation. This elevated precipitation results from orographic forcing, where moist air masses ascending the steep topography undergo adiabatic cooling, leading to enhanced condensation and rainfall. Other high-altitude stations, such as Mavrovo (<1,250 m), also exhibit high MAP values, reinforcing the role of elevation in amplifying precipitation through orographic effects.

Similarly, MAT values display significant spatial variability, primarily influenced by elevation and proximity to Mediterranean or continental climatic influences. The highest MAT, 15.5°C, is recorded in Gevgelija, a low-elevation station (<60 m) in the southeastern part of the country, which experiences a Mediterranean climate characterized by mild winters and warm summers. In contrast, the lowest MAT, 7.6 °C, is observed at Lazaropole and Mavrovo, both situated at high altitudes in the western mountains. This lower air temperature is attributed to the environmental lapse rate, approximately 6.5 °C per 1,000 m, which reduces air temperatures at higher elevations. Additionally, Lazaropole holds the record for the absolute minimum air temperature in North Macedonia, -31.2 °C, recorded during a winter period of intense radiative cooling under clear skies and stable atmospheric conditions. At higher altitudes, dynamic atmospheric circulation over steep slopes facilitates mixing and prevents prolonged stagnation of cold air, unlike in valleys and depressions, where cold air pooling can exacerbate low air temperatures.

As a major urban and industrial center, the climate of Skopje is significantly modified by the urban heat island (UHI) effect, a phenomenon where urban areas exhibit higher air temperatures than surrounding rural regions due to anthropogenic activities and infrastructure [76]. The UHI effect in Skopje, driven by heat-absorbing surfaces (e.g., concrete, asphalt), reduced vegetation, and waste heat from industrial and vehicular sources, elevates local air temperatures, particularly at night. Kaplan et al. (2018) [77] estimate that Skopje’s UHI can increase air temperatures by 1–3 °C compared to nearby rural areas, exacerbating aridity by enhancing evapotranspiration rates. This urban modification compounds the semi-dry conditions indicated by Skopje’s low I_DM, posing challenges for urban water management and heat stress mitigation.

Lazaropole, with an MAT of 7.6 °C, is the meteorological station with the lowest mean annual air temperature in North Macedonia, tied with Mavrovo. Its high-altitude location fosters colder conditions, and it is also the site where the absolute minimum air temperature of North Macedonia, -31.2 °C, was recorded. This extreme is attributed to radiative cooling in winter, amplified by clear skies and the station’s elevated, exposed position. In contrast to high-altitude stations, low-elevation valleys and depressions, such as those near Skopje and Štip, experience prolonged cold air pooling due to topographic sheltering, which can lead to air temperature inversions and locally lower minimum air temperatures. However, the dynamic circulation processes over steep slopes at higher altitudes, such as in Lazaropole, promote air mixing, reducing the persistence of extreme cold compared to valley bottoms. Thus, the western mountainous stations, exemplified by Lazaropole and Mavrovo, exhibit markedly more humid climates due to high precipitation and lower air temperatures, as reflected in their elevated I_DM values. In contrast, stations in the central and eastern lowlands, such as Štip and Skopje, are characterized by greater aridity, driven by lower precipitation, higher air temperatures, and, in Skopje’s case, urban-induced warming. These patterns align with regional topographic and climatic gradients, emphasizing the critical roles of elevation, continentality, and urbanization in shaping local aridity conditions across North Macedonia.

3.3. JCI distribution

The application of the Johansson Continentality Index (JCI) across North Macedonia reveals the presence of three of the four classical JCI climate classes—Moderate Marine, Transitional, and Continental—over both natural and agricultural areas (hereafter abbreviated as NatA and AgrA, respectively). The Marine class is notably absent, reflecting the minimal influence of maritime air masses within the national territory (Figure 5).

Moderate Marine climates are primarily observed in the western highlands and northeastern uplands, where topographic complexity and elevation act to moderate thermal extremes. Stations in this category include Mavrovo (30.1), Lazaropole (27.2), and Kriva Palanka (30.9). The reduced annual air temperature range in these areas reflects a substantial attenuation of continental effects, resulting in conditions partially influenced by distant maritime air masses (Table 6). Transitional climates occur predominantly in eastern peripheral and western regions, characterized by intermediate annual thermal amplitudes. Representative stations include Berovo (32.7) in the eastern highlands and Ohrid (32.4) in the southwestern region. These locations exhibit neither strongly continental nor distinctly marine characteristics, producing moderate seasonal air temperature variability. The Continental climate class is widespread, dominating the lowlands and central valleys of the country. Stations representing this regime include Bitola (39.4), Demir Kapija (39.7), Kumanovo (38.6), Skopje (38.8), Strumica (38.5), Štip (38.4), and Gevgelija (38.2). These areas are subject to pronounced seasonal air temperature contrasts, a consequence of their distance from maritime influences and the limited topographic shielding, which enhances the annual air temperature amplitude.

Analysis of the aggregated annual JCI revealed significant fluctuations in the thermal continental influence across the period of study. The highest JCI values were recorded in 2017 (approximately 49.33) and 2000 (approximately 47.05). These high index values indicate the most pronounced degree of continentality experienced in these specific years, reflecting maximum average annual air temperature amplitudes. Conversely, the lowest JCI value observed was 2014 (approximately 26.96). According to the categorization scheme adopted, this value places the year 2014 within the Moderate Marine climate category (JCI interval 20–32).

Thus, the JCI results demonstrate a distinct spatial gradient of continentality in North Macedonia. Moderate Marine climates are confined to the western highlands and northeastern uplands, Transitional climates are distributed in western peripheral and eastern border regions, while Continental climates dominate the central valleys and lowlands, encompassing the majority of the national territory. This spatial pattern emphasizes the integrative influence of latitude, topography, and distance from maritime sources in shaping the continentality regime across the country.

3.3. Monthly analysis

According to the monthly values of the I_DM (Im_DM), six out of seven climate types are recognized. January and December are considered extremely humid, with the maximal value of 71.9 (Figure 6). February, March, April, October, and November are very humid, where February has the highest value reaching 51.7. May is humid. September is characterized as semi-humid. June is Mediterranean, while July and August are semi-arid, both with a value of 14.6.

The I_DM reveals distinct spatial differentiation in climatic conditions across North Macedonia. The highest I_DM values, 72.9 at Lazaropole and 71.6 at Mavrovo, classify these locations as extremely humid (I_DM > 55), reflecting their high precipitation and cooler air temperatures driven by mountainous topography. These conditions support lush vegetation and water availability, characteristic of humid climates. In contrast, Štip (I_DM = 22.2) and Skopje (I_DM = 22.3) exhibit the lowest I_DM values, placing them within the Mediterranean climate category (20 ≤ I_DM < 24), indicative of semi-dry conditions. These lowland stations experience lower precipitation and higher air temperatures, contributing to greater aridity compared to the western highlands.

3.3. Seasonal analyses

The seasonal distribution of the Is_DM across North Macedonia reveals pronounced climatic variability controlled predominantly by topography, elevation, and regional meteorological influences (Figure 7).

Winter is generally humid across the country, with a national mean Is_DM of 58.7 (Table 7). Mountainous stations show extremely high humidity, with Mavrovo reaching 148.5 and Lazaropole 134.0, due to substantial precipitation (>300 mm) combined with low air temperatures (-1.2 to -0.7 °C). Lowland stations, including Skopje (35.4), Kumanovo (39.6), Prilep (34.3), and Strumica (46.4), exhibit semi-humid conditions, reflecting lower precipitation and warmer air temperatures (1–3 °C).

Spring displays generally humid conditions, national Is_DM 33.0, with the highest indices in mountainous stations: Mavrovo (66.8) and Lazaropole (65.1), supported by moderate precipitation (<270 mm) and rising air temperatures (<6–7 °C). Lowland areas such as Skopje (26.7), Kumanovo (26.6), and Strumica (26.0) experience semi-humid climates due to slightly reduced rainfall and warmer air temperatures (<12–13 °C).

Summer shows the highest climatic heterogeneity, national Is_DM 16.1. Mountainous stations maintain humid to semi-humid conditions: Mavrovo (20.9), Lazaropole (23.4), while lowlands experience pronounced aridity, for example, Gevgelija (12.3), Bitola (18.6), and Štip (17.9). High air temperatures (up to 27 °C) combined with reduced precipitation (<100–130 mm seasonally) explain the low summer indices in low-altitude stations. (see Table 8).

Autumn signals a return to humid conditions, national Is_DM 33.1. Mountainous regions such as Mavrovo (66.8), Lazaropole (70.8), Berovo (53.0), and Ohrid (49.6) show the highest seasonal indices. Lowland stations including Skopje (24.1), Kumanovo (25.9), Prilep (26.0), and Strumica (27.2) maintain semi-humid climates, reflecting moderate precipitation (<150–180 mm) and air temperatures (<12–14 °C).

This analysis emphasizes the critical influence of elevation, orography, and regional geography in shaping North Macedonia's climatic variability (see Figure 8). Humidity peaks consistently in northwestern mountainous regions during winter and spring, whereas lowlands and southern plains experience strong summer aridity. Autumn acts as a transitional season with spatial heterogeneity persisting. Such detailed seasonal patterns are crucial for hydrological heritage assessments, ecosystem management, agricultural planning, and climate adaptation strategies.

3.4. Annual De Martonne aridity index (I_DM) analysis

The De Martonne aridity index (I_DM) was calculated annually for 13 meteorological stations across North Macedonia, revealing significant spatial and temporal variability in hydrothermal conditions driven by altitude, precipitation patterns, and regional topography. Over the 30-year period from 1991 to 2020, annual I_DM values ranged from 18.8 to 42.3, reflecting a spectrum of climatic conditions from semi-arid to very humid, with a clear spatial pattern influenced by elevation and atmospheric dynamics. The mean annual I_DM across all stations, calculated as the average of the annual values, is 28.8, classifying North Macedonia’s overall climate as semi-humid on average, though with pronounced regional variations.

The highest I_DM values, indicative of the most humid conditions, are observed in the mountainous regions of western and northwestern North Macedonia. Mavrovo, with an I_DM of 75.7, and Lazaropole, with an I_DM of 73.3, both located in high-altitude areas near Mount Korab and the Bistra Mountain range, exhibit extremely humid conditions due to high precipitation and cooler air temperatures associated with mountainous terrain. Very humid climates dominate mid-altitude and mountainous regions, including Berovo with an I_DM of 44.6, Bitola at 42.1, Gevgelija at 37.1, Demir Kapija at 39.4, Kriva Palanka at 42.2, and Ohrid at 36.7, benefiting from consistent precipitation and moderate air temperatures (Table 9). In contrast, Strumica, with an I_DM of 29.4, exhibits a humid climate, reflecting lower precipitation and warmer air temperatures compared to higher-altitude regions. Semi-humid conditions prevail in lowland areas, including Kumanovo with an I_DM of 27.9, Prilep at 26.4, Skopje at 26.0, and Štip at 26.1, indicative of drier and warmer environments. Approximately 15% of the stations (two stations, Mavrovo and Lazaropole) are classified as extremely humid, 46% (six stations) as very humid, 8% (one station, Strumica) as humid, and 31% (four stations) as semi-humid, underscoring the dominance of humid conditions in North Macedonia’s climate.

The annual variations of the I_DM over the 30-year period show significant temporal variability (Figure 9). Of the 30 years analyzed, five years were classified as very humid (I_DM 35–55): 1995 with an I_DM of 36.5, 2002 with 40.3, 2009 with 37.2, 2010 with 42.3, and 2014 with 40.2. Sixteen years were classified as humid (I_DM 28–35): 1991 with 30.4, 1996 with 34.2, 1998 with 32.9, 1999 with 32.6, 2003 with 31.9, 2004 with 34.6, 2005 with 33.7, 2006 with 32.2, 2007 with 29.5, 2012 with 32.7, 2013 with 31.5, 2015 with 34.6, 2016 with 33.8, 2017 with 29.5, 2018 with 30.7 and 2020 with 28.4. Three years were classified as semi-humid (I_DM 24–28): 1997 with 26.8, 2008 with 25.2, and 2019 with 27.4. Five years were classified as mediterranean (I_DM (20–24): 1992 with 23.3, 1993 with 20.8, 1994 with 22.8, 2001 with 23.4, and 2011 with 20.8. One year, 2000, with an I_DM of 18.8, was classified as semi-arid, reflecting notably dry conditions consistent with regional drought patterns in the Balkans, as reported in studies such as Lukovic et al. (2013) [78] and Milevski et al. (2015) [79], and corroborated by the National Hydrometeorological Service – Skopje (2021) [80]. This dry year likely resulted from significantly reduced precipitation and elevated air temperatures. The maximum I_DM value was observed in 2014, with an I_DM of 42.3, aligning with reports of exceptionally high precipitation across the region, as noted by the National Hydrometeorological Service – Skopje (2021) [80] and the Copernicus Climate Change Service (2022) [81]. No years were classified as arid (I_DM < 10) or extremely humid (I_DM ≥ 55), indicating that North Macedonia’s climate remains predominantly within the humid to semi-humid range.

Although this analysis provides a robust understanding of the spatial and temporal variability of the I_DM, precise quantification of the area coverage for each climate type would require the development of a continuous annual I_DM map using spatial interpolation methods, such as Inverse Distance Weighting or Kriging. Such mapping would enable a detailed assessment of hydrological and ecological variability, supporting improved water resource management, climate impact analysis, and sustainable land-use planning.

Overall, the annual I_DM demonstrates a clear altitudinal control on climate, with elevated and mountainous regions maintaining extremely humid to very humid conditions, mid-altitude stations exhibiting very humid climates, and lowland areas displaying semi-humid characteristics (see Figure 10). The mean annual I_DM of 28.8, classifying the country as semi-humid, reflects the balance between the wetter, high-altitude northwest and the drier, lowland southeast. The rarity of semi-arid conditions, occurring in only one year, underscores the mitigating influence of North Macedonia’s continental and mountainous climate, which ensures a degree of moisture even in the driest periods. This pattern emphasizes the significant influence of topography and precipitation on the spatial distribution of climatic conditions in North Macedonia, providing critical insights for hydrological studies, ecological assessments, and climate adaptation strategies.

3.5. Annual Pinna combinative index (I_P) analysis

The Pinna combinative aridity index (I_P) was calculated annually for 13 meteorological stations across North Macedonia, revealing a diverse range of hydrothermal conditions shaped by topography and regional climatic factors. Over the 30-year period from 1991 to 2020, I_P values ranged from 9.4 to 37.7, reflecting a spectrum of climatic conditions from dry to humid, with a pronounced spatial pattern driven by elevation and atmospheric influences. The mean annual I_P for the entire country, calculated as the average of the annual I_P values over this period, is 18.2, indicating that North Macedonia’s overall climate is semi-dry (Mediterranean) on average, though significant spatial and temporal variations exist.

The highest I_P values, which indicate the most humid conditions, are observed in the mountainous regions of western and northwestern North Macedonia. Lazaropole, located in the high-altitude Bistra Mountain range, recorded the highest station-specific I_P at 39.5, followed closely by Berovo at 39.4 in the Maleševo Mountains and Mavrovo at 38.6 near Mount Korab, North Macedonia’s highest peak at 2,764 m. These stations experience substantial annual precipitation, such as 1058.5 mm at Lazaropole and 1046.0 mm at Mavrovo, combined with cooler mean annual air temperatures of 7.6°C at both locations, resulting in their classification as humid (I_P ≥ 20). The high elevation in these areas enhances orographic precipitation, as moist westerly air masses are forced upward, cooling and condensing to produce significant rainfall. The most humid year in the dataset was 2006, with an I_P of 37.7, likely most pronounced at Lazaropole due to its consistently high I_P values driven by elevation and precipitation.

In contrast, the lowest I_P values are found in the southeastern lowlands, where Demir Kapija recorded an I_P of 17.0, Gevgelija 18.2, and Strumica 19.1, all classified as semi-dry (Mediterranean) with I_P values between 10 and 20. These stations, situated in the Vardar River Valley near the Aegean Sea, experience higher mean annual air temperatures, such as 15.5°C at Gevgelija and 14.3°C at Demir Kapija, and lower precipitation in their driest months, for example, 25.8 mm in August at Gevgelija and 30.8 mm at Demir Kapija. The driest year recorded was 2000, with an I_P of 9.4, likely most pronounced at Demir Kapija, which consistently exhibits the lowest station-specific I_P due to its warm, dry conditions. This year marked the only instance of fully dry conditions (I_P < 10) in the 30-year period, likely driven by anomalously low precipitation and elevated air temperatures influenced by strong Mediterranean air masses. The lowest I_P value, 9.4 in 2000 at Demir Kapija, reflects the extreme aridity in the southeastern lowlands during that year.

Central and northern stations, including Kumanovo with an I_P of 25.4, Skopje at 23.9, Ohrid at 21.5, and Štip at 21.5, are classified as humid, with I_P values ranging from 21.5 to 25.4. Prilep, with an I_P of 30.6, and Kriva Palanka, with an I_P of 35.7, also exhibit humid conditions, the latter driven by substantial annual precipitation of 661.1 mm despite a moderate mean air temperature of 10.7°C. Bitola, with an I_P of 20.3, lies in a transitional zone between the humid northwest and semi-dry southeast.

Over the 30-year period, climatic conditions varied significantly (Figure 11). Dry conditions, defined as I_P values below 10, occurred only in 2000 with an I_P of 9.4, representing just one year or 3.3% of the period. This exceptionally dry year, likely most pronounced at Demir Kapija, was driven by significantly reduced precipitation and higher air temperatures, possibly due to strong Mediterranean air masses affecting the southeastern lowlands. Semi-dry conditions, with I_P values between 10 and 20, were observed in 20 years, accounting for 66.7% of the period: 1992 with an I_P of 11.7, 1993 with 10.4, 1994 with 13.4, 1995 with 18.2, 1996 with 18.3, 1997 with 14.2, 1999 with 18.0, 2001 with 13.0, 2003 with 16.3, 2007 with 14.7, 2008 with 13.1, 2011 with 10.5, 2012 with 16.3, 2013 with 16.4, 2015 with 17.2, 2016 with 16.9, 2017 with 17.8, 2018 with 15.3, 2019 with 14.0, and 2020 with 14.5. These years were characterized by warmer, drier Mediterranean influences, particularly affecting southeastern stations like Demir Kapija, Gevgelija, and Strumica. Humid conditions, with I_P values of 20 or greater, prevailed in 9 years, representing 30.0% of the period: 1991 with an I_P of 25.5, 1998 with 20.3, 2002 with 25.5, 2004 with 25.4, 2005 with 27.6, 2006 with 37.7, 2009 with 28.0, 2010 with 21.2, and 2014 with 21.1. These years corresponded to periods of higher precipitation and cooler air temperatures, particularly in the mountainous northwest, with stations like Lazaropole and Mavrovo experiencing the most pronounced humidity.

The spatial distribution of I_P values reveals a clear gradient from the humid northwest to the semi-dry southeast. High mountains in the west, such as Korab, Bistra, and Šar Planina, act as barriers that capture moisture from westerly air masses, resulting in higher precipitation and cooler air temperatures. Conversely, the southeastern lowlands, influenced by the Vardar valley’s channeling of warm, dry Mediterranean air from the Aegean Sea, experience reduced precipitation and higher air temperatures, leading to semi-dry conditions. The Aegean Sea’s proximity exacerbates aridity in the southeast, particularly in summer months, as evidenced by low precipitation values. The mean annual I_P of 18.2, classifying North Macedonia as semi-dry on average, reflects the balance between the humid, high-altitude northwest and the drier southeast. The rarity of fully dry conditions, occurring in only one year, underscores the mitigating influence of the country’s continental and mountainous climate, which ensures a degree of moisture even in the driest regions (Figure 12). The marked humidity at Lazaropole, Mavrovo, and Berovo emphasize the critical role of elevation, while the semi-dry conditions in the southeast, particularly at Demir Kapija, reflect the Aegean Sea’s influence, creating a distinct climatic gradient across North Macedonia (Table 10).

3.6. Aridity trend

Aridity trends in North Macedonia were analyzed using the JCI, I_DM, and I_P indices through three statistical methods. First, linear trend equations were derived from mean annual JCI, I_DM, and I_P values across the region, graphically illustrating aridity changes [43,82,83]. Second, the Mann-Kendall test assessed significant trends, testing the null hypothesis (H₀: no trend) against the alternative (H_a: significant trend) at α = 0.05. A p-value < α rejects H₀, confirming a trend, while p ≥ α retains H₀ [84,85,86,87,88]. This robust approach, effective for missing data [43], enhances understanding of North Macedonia’s climatic trends. Third, to refine the analysis and detect possible shifts in aridity over the study period, Pettitt’s nonparametric test was employed.

The temporal variation of aridity in North Macedonia was further analyzed using linear regression. For the I_DM, the obtained regression equation was: y = 0.028x − 25.65 and for the I_P: y = −0.011x + 37.15, where t is the time in years. For the JCI, the obtained regression equation was: y = 0.1514x +33.78. Trend lines are shown in Figure 13. I_DM equation indicates a weak positive tendency in the I_DM, while I_P equation reflects a very weak negative tendency in the I_P, and JCI equation shows a very weak positive tendency in continentality. In all cases, the trends are extremely weak.

The Mann–Kendall (MK) test was applied to determine the statistical significance of the observed tendencies (Table 11). As the obtained probabilities were greater than the significance level (α = 0.05, i.e. 5%), the null hypothesis H₀ of no trend cannot be rejected. The risks of committing a Type I error (rejecting H₀ while it is true) were 48.6% for the I_DM, 85.8% for the I_P, and 72.4% for the JCI, respectively. Consequently, the MK test showed that, considering all values of both indices, there is no statistically significant change in aridity trends in the last 30 years in North Macedonia.

The table summarizes Mann-Kendall trend analyses for the JCI, I_DM, and I_P, across 13 meteorological stations in North Macedonia (1991–2020), reporting annual averages and Sen's slope (β) estimates of long-term change. It shows no significant trends in aridity or continentality (all Z_MK < 1.96, p > 0.05), underscoring overall climatic stability, with minor negative β values in low-elevation sites indicating subtle drying tendencies. This spatial pattern—characterized by wetter mountainous areas and drier, more vulnerable valleys—provides an essential reference point for understanding how ongoing climate warming may influence agricultural productivity and water availability across the region.

3.7. Correlation

Analysis of the I_DM and I_P for North Macedonia from 1991 to 2020 reveals stable climatic conditions, with no significant trend in aridity, indicating consistent precipitation and air temperature patterns across the region (Figure 14). Both indices demonstrated comparable values in spatial distribution and no notable change in aridity, prompting a correlation analysis between them. The Pearson correlation coefficient r has a value of 0.66. The coefficient of determination r2 is 0.437 when the regression type is linear, as shown in Figure 14, with the equation y = 0.691x – 3.108.

Correlation coefficients between the North Atlantic Oscillation Index – NAO [89], El Niño-Southern Oscillation Index – SOI [90], and seasonal and annual I_DM were determined to investigate possible relationships between aridity and atmospheric variability. To ensure transparency and reproducibility, NAOI and SOI data were sourced from NOAA's Climate Prediction Center (for NAOI) (National Centers for Environmental Information, 2023) and Physical Sciences Laboratory (for SOI) [91], with monthly values averaged to derive seasonal indices (Winter: December–February, Spring: March–May, Summer: June–August, Autumn: September–November).

Annual correlations were derived by averaging seasonal I_DM, NAOI, and SOI values for each year. Results are shown in Table 12. The correlation between the I_DM and the NAOI is significant (at the 5.0% level of significance) in summer (0.45). The positive correlation coefficient in summer indicates a drying effect on the I_DM in North Macedonia, potentially due to reduced precipitation under positive NAOI phases. A negative correlation coefficient in winter (−0.31, p = 0.098) suggests a wetting effect on the I_DM, though not statistically significant, reflecting increased moisture during negative NAOI phases. Since precipitation is one of the key components of the I_DM, the response of the I_DM to the NAOI can be associated with the established links between this index and precipitation dynamics [92]. Our results align with findings from regional studies on summer precipitation patterns in the Southeastern Europe [93]. For the SOI, a weak positive correlation in autumn (0.19) suggests minor ENSO influence during spring-like transitions in North Macedonia.

Coefficients significant at the 5% level (α = 0.05), indicating less than a 5% probability of Type I error, are presented in bold

Finding similar or identical results at regional and national scales is challenging and difficult. However, certain spatio-temporal patterns of aridity can be identified. This difficulty largely stems from the nature of the climate indices chosen for the study, as well as the specific selection criteria employed by the authors.

The spatial and temporal patterns of continentality and aridity in North Macedonia, derived from the JCI, I_DM, and I_P, reveal semi-dry conditions in lowland southeastern areas and semi-humid conditions in highlands during summer. These patterns are consistent with observations in neighboring countries. In Bulgaria, analysis of the I_DM index during 1961–2015 revealed semi-arid to arid conditions, while multiannual variability indicated decreasing aridity over the last three decades [53]. Serbia shows pronounced continentality in lowlands with semi-dry to semi-humid climates and an absence of Mediterranean-type dry zones [45]. In Kosovo and Metohija, semi-dry and semi-humid climates dominated during 1965–1999, with no significant aridity change [44]. Montenegro exhibits generally humid to extremely humid conditions, except for semi-arid summer periods in coastal and adjacent inland areas, with trends indicating increasing aridity in April and August over the last sixty years [47]. North Macedonia’s lowland summer aridity and highland semi-humid conditions correspond closely to these regional observations, supporting the consistency of hydroclimatic patterns across the region.

Expanding to South-Eastern Europe, Turkey exhibits semi-arid and dry sub-humid conditions in continental interiors and Mediterranean regions, with decreasing annual and winter precipitation since the 1970s, and seasonal increases in spring precipitation in some areas [34]. Iraq, in contrast, is characterized by arid and semi-arid conditions covering 97% of the country, with increased dryness after 1997, particularly in northern regions [40]. These findings reflect the sensitivity of lowland North Macedonia to summer water deficits, highlighting comparable responses to hydroclimatic variability in continental interiors and semi-arid zones.

Projections for South-Eastern Europe (1961–2050) based on multiple climate indices, including the JCI, I_DM, and I_P, indicate that the overall spatial pattern of aridity remains relatively stable, yet marked shifts toward more arid conditions are expected in the Balkan Peninsula, including North Macedonia [65]. These large-scale trends are corroborated by ultrahigh-resolution analyses of regional bioclimate. Charalampopoulos et al. (2023) [57], using the I_DM, projected future bioclimate changes across natural and agricultural areas, highlighting a pronounced expansion of semi-dry conditions in central, southern, and southeastern parts of North Macedonia under high-emission scenarios for the period 2071–2100. In addition, Charalampopoulos and Droulia (2025) [58] reported increasing aridity and continentality patterns over South-Eastern Europe from 1964–2023, with semi-dry conditions expanding notably over agricultural and natural areas, particularly in Bulgaria, Serbia, and North Macedonia. The results of the present study align with these projections, demonstrating that lowland areas in North Macedonia are increasingly affected by summer aridity, whereas highlands remain comparatively buffered. This agreement underscores the robustness of the applied indices for capturing hydroclimatic gradients and supports the broader regional relevance of observed seasonal and spatial aridity trends.

Regional studies across the South-Eastern European and Mediterranean domains indicate that topography, latitude, and continentality strongly modulate seasonal aridity. In Southern Bulgaria, decreasing multiannual aridity trends were observed [53], whereas Montenegro’s coastal summer zones exhibited semi-arid conditions [47]. Southwestern Spain experienced increasing aridity between 1951 and 2010, particularly in historically humid northern regions, with spring and summer showing the highest water deficits [94]. Italy shows a progressive reduction in precipitation and increasing temperatures between 1951 and 2010, leading to a widening gap between water demand and rainfall in agricultural systems [30]. North Macedonia similarly displays summer lowland aridity, while highlands retain semi-humid conditions, confirming the generality of seasonal patterns across Mediterranean and South-Eastern European terrains. Furthermore, Greece demonstrates a gradual shift from moderate marine climates to more continental and semi-dry classes in northern and central regions, while islands remain stable [95]. Comparable continentality–oceanity gradients were observed in the Iberian Peninsula, with maritime northern areas and continental interiors, reflecting similar topographic and latitudinal modulation of climate as seen in North Macedonia [96]. These observations corroborate the influence of topography and latitude on spatial aridity variability, which aligns with the patterns observed in North Macedonia.

Moving eastward, Central Kazakhstan exhibits a northeast–southwest aridity gradient, with southern desert areas becoming drier and northern highlands more humid between 1960 and 2022 [97]. Seasonal analysis revealed winter moisture gains contrasted by extreme summer aridity, consistent with North Macedonia, where lowlands are semi-dry in summer, and highlands maintain semi-humid conditions. Similarly, recent analyses across the Middle East and North Africa [98] indicate increasing continental influences and modest aridity trends, particularly over humid regions such as the Ethiopian Highlands and Fertile Crescent, underscoring the role of large-scale circulation patterns in modulating spatial hydroclimatic.

At the European scale, continental and Mediterranean regions display heterogeneous aridity trends. In Italy, high spatial variability in aridity and continentality indices reflects mountain-induced precipitation gradients [99]. Across Europe, long-term datasets indicate expansion of dry regions since the 1980s, particularly at low and medium altitudes in Western, Central, and Eastern Europe [100]. According to high-resolution projections of aridity under climate change, the South Balkan Peninsula, including North Macedonia, experiences dry climates across all indices, with southern lowlands classified as semi-dry and highlands maintaining semi-humid conditions. Our findings in this study overlap with these continental-scale projections, confirming the reliability and regional representativeness of the applied indices.

Yet, certain limitations are inherent in this study. The analysis is restricted to the period 1991–2020, which constrains the assessment of longer-term trends and decadal variability. While station-based observations provide high-resolution and reliable data for the examined period, extending the temporal coverage to include earlier records or climate projections would allow a more comprehensive evaluation of hydroclimatic changes. One of the limitations is the small number of meteorological stations (13) used in the study, leading to a high degree of generalization and interpolation. Furthermore, the integration of satellite-derived datasets, such as MERRA-2 or ERA5 reanalysis products, alongside geospatial platforms like Google Earth Engine (GEE), would enable cross-validation, spatial extrapolation, and the construction of broader temporal and spatial scenarios [4]. Incorporating these approaches would enhance the robustness and representativeness of regional aridity assessments, thereby supporting improved water resource management and climate adaptation strategies [5,101].

Thus, the North Macedonian patterns of semi-dry lowlands and semi-humid highlands, with summer as the critical period for water deficits, are broadly consistent with regional and extra-regional studies across South-Eastern Europe, the Mediterranean, and Central Europe and Asia. This congruence validates the applicability of the JCI, I_DM, and I_P, indices for regional hydroclimatic assessment and underscores the relevance of these findings for water resource management, agriculture, and climate adaptation strategies.

4. Conclusions

In this study, the long-term spatio-temporal variability of aridity indices in North Macedonia from 1991 to 2020 was analyzed. For the purpose of the study datasets from 13 meteorological stations were used. Certain patterns can be identified so general conclusions can be summarized as follows:

Mean annual temperatures and mean annual precipitation sums showed spatial variability influenced by hypsometry and the proximity of major water aquatories such as Mediterranean basin. Therefore, the study area has diverse climate (from continental to mild, Mediterranean influences). However, the urban agglomeration of Skopje is affected by UHI phenomenon which leads to increasing of mean annual temperature. This fact is the cause of pronounced continentality and semi-dry conditions.

On a temporal scale the highest JCI values were recorded in 2017 and 2000 which suggest high continentality. The lowest JCI value observed was during 2014 (Moderate Marine climate class). The spatial distribution of the JCI detected three classes – moderate Marine, Transitional, and Continental climate. Moderate Marine climates are observed in the western highlands and northeastern uplands. Transitional climates is present in eastern and western regions while continental climate is dominant in the central parts of the country and in central valleys.

The monthly I_DM values show the increasing of aridity from June to September. July and August belong to semi-arid classes while June is Mediterranean. This fact could have influence on agriculture in North Macedonia and create a need for expanding areas under irrigation. The seasonal values of this index indicate extremely humid climate during winter season, spring and autumn are characterized by a humid climate, while summer is characterized as semi-arid. In term of spatial changes, summer season indicates semi-arid conditions in most of North Macedonia. The annual values of the I_DM indicate the following: 5 years were very humid, 16 humid, 3 semi-humid, 5 mediterranean and only one year (2000) belongs to the semi-arid category. In general, according to the I_DM, North Macedonia has humid to semi-humid conditions on annual scale. Spatial distribution of the I_DM indicates dominant humid, very humid and extremely humid classes.

Mean annual value of the I_p reveals dominant semi-dry (Mediterranean) influences – 20 years. Humid classes are pronounced during nine years. Similarly, just one year based on the Pinna combinative index belongs to the semi-dry conditions (2000). Spatial distribution of I_p detected humid class on almost whole meteorological stations in the country. However, three meteorological stations (Gevgelija, Demir Kapija and Strumica) showed different aridity patterns – semi-dry (Mediterranean) climate.

The MK test indicates no pronounced statistical significance in the aridity trends. The Pearson correlation test between the I_DM and the NAOI is statistically significant in summer. The positive correlation in summer indicates a drying effect on the I_DM in study area (linked to reduced precipitation under positive NAOI phases). A negative correlation coefficient in winter displays a wetting effect on the I_DM (no statistical significance) reflecting increased moisture during negative NAOI phases.

The results of this research provide a solid basis for future studies of aridity in North Macedonia and can be of importance in numerous aspects of anthropogenic activities. Here, policy planning on a national scale, sustainable agricultural production and sustainable management of water resources (monitoring of irrigation conditions) are emphasized. Considering that increasing of aridity can be expected in the future, this will have a negative impact on traditional agriculture, which is the most widespread in the study area. Consequently, state and local authorities need to implement measures and develop plans for sustainable agriculture, particularly for crops reliant on summer precipitation (e.g. using drought-resistant crops).

Future investigations of aridity variability in North Macedonia must include a wider set of different aridity indices or implementation of remote sensing techniques in aridity detection. Also, a long-term assessment of aridity changes is desirable (by expanding used time series to the historical data). Since the aridity conditions are related to drought conditions, this phenomenon should be investigated by using the appropriate empirical methods such as: for meteorological drought (Standardized Precipitation Index, SPI), Standardized Precipitation-Evapotranspiration Index, SPEI), for hydrological drought (Palmer Hydrological Drought Index, PHDI, Streamflow Drought Index SDI), for agricultural drought (Crop Moisture Index, CMI and Vegetation Condition Index, VCI). This approach provides the opportunity of comparing results between the aridity and drought and could give a valuable insight into the moisture deficit and its dynamics on a national scale.