1. Introduction

Climate change and climate variability represent a risk for agriculture, food security, and the Colombian economy and society (^{Hodson et al., 2017}). Evidence of climate change in Colombia has been described in several studies such as ^{Pabón-Caicedo (2003)}, ^{Vuille et al. (2003)}, ^{Ochoa and Poveda (2008)}, ^{Cantor-Gómez (2011)}, ^{Carmona and Poveda (2014)}, and ^{Hurtado-Montoya and Mesa-Sánchez (2015)}. ^{Ramírez-Villegas et al. (2012)} warn that by 2050 climate change will affect about 3.5 million people in Colombia, affecting 14% of the gross domestic product (GDP), corresponding to agriculture, employment, agro-industries, supply chains, and food and nutrition security. It is a priority to understand annual and interannual climate variability and its correlation with weather events in order to mitigate their negative impacts on national agricultural production (^{Barrios-Pérez et al., 2021}).

Climate behavior based on meteorological variables such as precipitation and air temperature has been previously studied (^{Pabón-Caicedo et al., 2001}). Studies in Colombia show that interannual variability of temperature responds to factors such as El Niño-Southern Oscillation (ENSO) (^{Poveda, 2004}). ^{Puertas-Orozco and Carvajal-Escobar (2008)} evidence the correlation between air temperature and precipitation during ENSO events, with precipitation decreasing as air temperature increases, and also show that temperature varies with the physiography of the country.

In this sense, temperature correlates with NIÑO 3.4 and NIÑO 4 during December-January-February, and these correlations improve when considering the interannual variability (^{Puertas-Orozco and Carvajal-Escobar, 2008}). ^{León (2000)} shows positive trends in maximum recorded temperatures, decreasing trends in mean temperature, and evidence of spectral frequencies with ENSO events, but this study used just one meteorological station in the northern Colombian region.

Mean annual temperatures in the Colombian Caribbean region range between 28 and 30 ºC, except in the Sierra Nevada de Santa Marta (SNSM). The annual temperature cycle allows us to identify three seasons: (1) the warm season, with maximum temperatures from February to April, (2) the transition season, from June to August, and (3) the cold season, with the lowest temperatures from September to November (^{Pabón-Caicedo et al., 2001}).

Considering the evidence of climate change in Colombia, its impacts on agriculture, the economy, the health sector, food security, and the fact that studies carried out in northern Colombia used few stations to describe the spatial and temporal patterns of air temperature, it is necessary to evaluate the variability, cycles, and trends of mean air temperature in this region. Therefore, using historical temperature series from 24 meteorological stations distributed in three departments located in northern Colombia, this study seeks to contribute to the understanding of air temperature variability patterns in the Colombian Caribbean, the periods of oscillation of this variable, its relationship with climate variability indices, and temporal trends for a period of 30 consecutive years until 2020.

2. Study area

The Colombian Caribbean is divided into eight regions or departments (^{Meisel-Roca and Pérez-Valbuena, 2006}). Seven of them are located in the continental Caribbean plain ending in the SNSM, which descends in the Guajira peninsula (^{Rangel-Buitrago and Anfuso-Melfi, 2013}). The SNSM joins three departments in the northeastern Colombian Caribbean: Magdalena, Cesar, and La Guajira (between 8-11º N and 72-74º W) (Fig. 1), giving it distinctive features of relief and climate compared to the rest of the region (Meisel-Roca and Pérez-Valbuena, 2006). According to the Caldas-Lang classification of ^{IDEAM (2015)}, the northern Colombian Caribbean is characterized by a dry climate in the low and flat parts. In mountainous areas, the climate varies from warm and humid to temperate and cold (^{IDEAM, 2011}). In La Guajira, the northernmost part of South America, the climate is arid (^{IDEAM, 2011}; ^{Vanegas-Chamorro et al., 2015}).

Fig. 1 Study area. Departments of La Guajira (purple), Cesar (yellow), and Magdalena (orange), and meteorological stations (red triangles). The description of the stations can be found in Table I. 

In the Department of Magdalena, the northeastern trade winds interact with the SNSM, creating a valley-mountain circulation that causes temperatures to drop in December, January, and February. In March, the heat is accentuated and humidity increases. In La Guajira maximum temperatures exceed 34 ºC in June and August (^{Martín-Garzón, 2021}). In the department of Cesar the average annual temperature exceeds 28 ºC. Towards the foothills of the SNSM and the Serranía del Perijá (SP), the temperature has a relationship with altitude, which characterizes the Colombian mountains (^{IDEAM, 2015}). The rainfall regime in Magdalena and Cesar is bimodal, with the first rainy season from April to May and the second one, more intense, from September to November. In La Guajira the rainfall is scarce and unevenly distributed, oscillating around 300 mm yr^-1 (^{IDEAM, 2015}).

3. Methods

Monthly mean air temperature data (ºC) from 24 meteorological stations (Fig. 1) were provided by the Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM), available on the website http://dhime.ideam.gov.co/atencionciudadano/. Table I describes in detail the characteristics of the stations used.

A quality control was carried out to verify the reliability of the datasets, considering the number of missing data in each of the time series and their temporal length. Time series with three or more months missing were excluded. As a result, 24 were available for the departments of Magdalena, Cesar, and La Guajira. Nonlinear principal component analysis (NLPCA) was used to complete the missing data in each time series, parametrized with a circular inverse model. The network is trained in 1500 iterations (^{Scholz et al., 2005}).

A circular unit at the component layer is introduced to the NLPCA, which describes a potential circular data structure by means of closed curve, to obtain circular components. The auto-associative neural network contains a circular unit pair (p, q) in the component layer, and the output values z _p and z _q (Eq. [1]) are constrained to lie on a unit circle represented by a single angular variable θ (^{Scholz, 2007}).

The forward propagation is represented by:

where z _m are sums of their inputs given weights (ω _pm and ω _qm ) for the matrix. For backward propagation, we need the derivatives of the error function remaining for all n samples to be:

With the reconstructed series, the annual and semi-annual variability was analyzed using monthly climatology and boxplot diagrams.

To evaluate the performance of the missing data reconstruction method (NLPCA), we used the Spearman correlation coefficient (r_sp) (Eq. [4]) (^{Spearman, 1904}; ^{Zar, 1972}) with a 95% significance. Mean absolute error (MAE) (Eq. [5]), and bias (Eq. [6]) are employed in the paired comparisons between the reconstructed and original time series.

In these equations, D is the difference between the ranks of two samples, N the number of ranks, ts _r is the reconstructed time series, ts _o represents the original time series, and n is the total number of data.

The critical value (r _c ) for the Spearmean correlation (Eq. [7]) is calculated from the following expression:

where the hypothesis test is given by:

The WMO recommends using the most recent 30-year period ending in a year ending in 0 for the calculation of the climate normal for the current period from 01/01/1991 to 31/12/2020 (^{WMO, 2017}). Therefore, the time series are cut for this period. The time series periodicities were evaluated using the cyclic descent technique (^{Bloomfield, 1996}), and implemented with the tool developed by ^{González-Rodríguez et al. (2015)}. This technique consists of decomposing the time series to detect the most important frequency, phase, and amplitude, as well as statistically significant harmonics.

where α and β are the linear regression coefficients when a linear trend is assumed (otherwise both parameters equal zero), and a _i , b _i and p _i are the corresponding parameters of the i _th harmonic. Thus, a _i and b _i are used to calculate the amplitudes and phase angles.

If a linear trend exists in the time series, its magnitude can be quantified using ^{Sen’s (1968)} non-parametric slope estimator:

where β represents the slope mean values between data measurements x _j and x _i at time steps j and i, respectively, with j > i. The positive value of β indicates an increasing trend while its negative value indicates a decreasing trend. The direction of the data trend is indicated by the sign of β, while its slope is determined by its value. By comparison to linear regression, this technique offers the benefit of constraining the impact of missing data or outliers on the slope.

To establish a relationship between the mean temperature at the 24 selected stations and some climate variability indices, non-parametric correlations were performed using Spearman’s coefficient (r _sp ) (Eq. 4) and the hypothesis test given in Eq. (8). The variability indices used are the Southern Oscillation Index (SOI), which is a standardized index based on the observed sea level pressure (SLP) differences between Tahiti and Darwin, Australia (^{NOAA, 2005a}); the Multivariate ENSO Index (MEI v2), which is a bi-monthly time series that represents the combined empirical orthogonal function (EOF) of five variables: SLP, sea surface temperature (SST), zonal and meridional components of surface wind, and outgoing longwave radiation (OLR) over the tropical Pacific basin (30º S-30º N, 100º E-70º W) (^{Wolter and Timlin, 2011}; ^{NOAA, 2018}). The North Atlantic Oscillation (NAO) index relies on the surface SLP difference between the Subpolar Low and the Subtropical (Azores) High. A positive NAO phase signifies pressure and height anomalies that are below average across the North Atlantic high latitudes and above average across the central North Atlantic, western Europe, and the eastern United States (^{NOAA, 2005b}). The sunspot number is an index of the activity of the entire visible disk of the sun. It is determined each day without reference to preceding days and determined by the sunspot group, and it may consist of one or a large number of distinct spots of different sizes (^{Künzel, 1962}; ^{SILSO, 2015}).

Two non-parametric techniques, the modified Mann-Kendall test (^{Hamed and Rao, 1998}) and the Theil-Sen estimator (^{Theil, 1950}; ^{Sen, 1968}), also used in works such as ^{Dahmen and Hall (1990)} and ^{Morales-Acuña et al. (2019)}, were used to calculate the monthly mean temperature trend and multi-year monthly averages.

The modified Mann-Kendall test offers an alternative for detecting trends in time series that have autocorrelation with the difference in the formula to calculate variance (^{Morales-Acuña et al., 2021}):

where Var(S)* is the modified variance, n the number of observations, m the number of data pairs in the series, and t _i the number of ties at the timestep i. Thus n/n* (Eq. 12) is the correction factor given by r _k ^R , corresponding to the lag-k autocorrelation coefficient of the ranks of the data:

The significance of the trend is determined with Eq. (13):

4. Results

The mean air temperature time series were reconstructed with 95% reliability to complete the missing data, with a mean correlation coefficient of 0.99 and a mean absolute error of 1.8 × 10^-4 ºC. The performance of this technique is shown in Figure 2.

Fig. 2 Reconstruction of time series by the NLPCA method and statistical metrics between the original and the reconstructed time series of RINCON station (department of Cesar). The results of the time series reconstruction are shown in full in the supplementary material (Figs. S1, S2 andS3). 

The analysis of the annual cycle (Fig. 3) allowed us to identify that the stations located in the department of Cesar (Fig. 3a-i) are characterized by a bimodal behavior, in which the stations have their highest maximum temperature in March and the second maximum temperature occurs in July. Temperature ranges from 27 to 31 ºC. Likewise, the minimum of maximum temperatures is found at the RINCON station (350 masl) (Fig. 3h) located at the foothills of the SP.

Fig. 3 Monthly climatology of air temperature in the departments of (a-i) Cesar, (j-r) La Guajira, and (s-y) Magdalena for the period 1991-2020. The green line is the monthly mean. Red crosses are outliers. 

The behavior in La Guajira is mostly unimodal (Fig. 3j-r). The highest temperatures occur mostly in June-July-August; however, the department can be divided into two zones: north and center-south. In the first one, all stations obey this behavior except NAZARE (Fig. 3o), which has a minor first high temperature in June. Also, URUMIT (Fig. 3r) in the south of the department has a similar behavior to stations in the Cesar department, since it has two maximum temperatures in April and July (approximately 29 ºC).

The department of Magdalena (Fig. 3s-y) is divided into north, center-south, and mountain zones (Fig. 3y). The north-center fraction is formed by a station in the coastline (ASIMBO, Fig. 3s), which is unimodal with a maximum temperature of ~ 29 ºC in May. The station with highest altitude and lowest maximum temperature is SLOREN (2200 masl) (Fig. 3y), located in the foothills of the SNSM; it has a unimodal behavior with maximum temperature in July and minimum high temperature. In this station, the super humid cold weather predominates, and temperatures fluctuate between 12.5 and 14.5 ºC. The center-south fraction has a bimodal behavior, within which the stations located at the center (PADELM, Fig. 3v, and MELUNA, Fig. 3u) have their first highest maximum in April and the second in July. The stations further south in the department (SEISEL, Fig. 3x, and LALAMO, Fig. 3t) have a highest maximum temperature in March of ~ 29 ºC, and a second maximum temperature in July with an average of 27.4 ºC.

The periodicity analysis (Table II) allowed us to identify from the five significant harmonics that air temperature variability responds mainly to the annual cycle, except for CHIRIG, located in the valley of Cesar, where the highest significance corresponds to a ~ 6-years cycle and the annual cycle in the fifth harmonic. For the second harmonic, the stations are influenced by semi-annual periods (22%), and periods of 3-7 years (16%) and more than 12 years (28%). In the third, fourth, and fifth harmonics, 3-7 years periodicities predominate with 38%, 56%, and 34% of occurrence, respectively.

For all stations, the Spearman correlation mean coefficient with 95% of significance and N = 360 is greater than the critical value (r _c = 0.103) rejecting the null hypothesis. In this respect, it is shown that there is a better association with the MEI (83%), SOI (75%), and sunspot (25%).

The correlation with the MEI (Fig. 4b) is mostly significantly moderate in Cesar (13%) and La Guajira (13%), with a correlation coefficient of r_sp = 0.66. The northern region of Cesar and central Magdalena have a high positive correlation (4%) and a moderate positive correlation (17%). The correlation with the SOI (Fig. 4a) is inverse, mostly with low correlation in Cesar (17%) and La Guajira (17%) with r_sp = -0.28 and moderate correlation in Cesar (13%) and La Guajira (13%), with r_sp = -0.39. The strongest signals of correlation with SOI are in the coast of Magdalena, SNSM, and the north part of Cesar.

Fig. 4 Spearman’s correlation coefficient (95% significance, and critical value of ± 0.103) between air temperature and (a) SOI, (b) MEI, (c) NAO, and (d) Sunspot index. 

It can be said that with 95% confidence there is no relationship with the NAO (Fig. 4c). The negative correlation that exists is very low, localized in the coast of Magdalena, in the northern part of La Guajira, and near Urumita. A similar situation occurs with sunspots (Fig. 4d) and mean air temperature, since there is no relationship. The signal is negative very low and occurs in the foothills of the SNSM towards La Guajira.

The trend in all stations is positive, indicating an increase in temperature. The maximum trend in full series (Fig. 5) is recorded in the department of Cesar at the GUAYMA station, with 0.07 ºC yr^-1. This is followed by the maximum in the department of Magdalena at the PADELM station, with 0.05 ºC yr^-1, both between 20-50 masl. Lesser positive trends (0.02 ºC yr^-1) in ASIMBO and LALAMO (Magdalena), and 0.01 ºC yr^-1 in RANCHO (La Guajira) are found.

Fig. 5 Trend analysis (ºC yr^-1) of average temperature of the full series. 

Similarly, the monthly multiannual trend analysis (Fig. 6) shows positive values in all months, ranging between 0.2and 1.07 ºC yr^-1. The maximum positive trends occur most of the year in the center of Cesar (GUAYMA reaching 1.07 ºC yr^-1 in June); however, January (0.72 ºC yr^-1) and August (0.68 ºC yr^-1) the maximum positive trends occur in PADELM (northeast of Magdalena). The minimum trends occur in March in RANCHO (0.2 ºC yr^-1). Minimum trends occur in January near the coast of Magdalena (ASIMBO, 0.32 ºC yr^-1), and between February and May, July-August and October-November in mountain foothills. In the north coast of La Guajira, the minimum trends occur in the PBOLIV station (June and September, with 0.28 and 0.38 ºC yr^-1, respectively) and NAZARE (December, with 0.27 ºC yr^-1).

Fig. 6 Trend analysis (ºC yr^-1) of monthly average temperature. 

5. Discussion

The reconstruction of mean air temperature time series from the inverse NLPCA approach (Fig. 2) allowed us to identify that the extraction of the nonlinear components from the incomplete time series and its subsequent nonlinear mapping achieved an adequate reconstruction of the analyzed time series. Although the reconstructed series slightly overestimate the mean values, the variability of the original series is preserved, as reported by other authors when reconstructing time series of different variables (^{Ruessink et al., 2004}; ^{Barakzehi et al., 2013}; ^{Canchala et al., 2020}; ^{Morales-Acuña et al., 2021}).

The analysis of monthly air temperature climatologies in the Colombian Caribbean (Fig. 2) identified a heterogeneous distribution in terms of mean air temperature values, such that a warm season with maximum temperatures in July was identified in the department of La Guajira, while two warm seasons with maximum temperatures in April and July were found in the department of Cesar. In Magdalena, due to the presence of the SNSM, the proximity to the Caribbean Sea, the extensive plains located in the center of the department, and its interaction to the south with the La Mojana region, a warm season with maximums in May and July was characterized for the coastal region and the SNSM, while in the plains of the central region two warm seasons with maximums in April and July were observed and in the southern region two warm seasons with maximums in March and July were identified. This heterogeneous distribution of air temperature responds to regional differences due to latitude, topography, mountain range orientation, and the average vertical temperature gradient (5.53 ºC km^-1) reported by ^{Eslava (1992)} for all of Colombia. In addition, the seasons and heterogeneity found in the present study coincide with those reported by ^{Pabón-Caicedo et al. (2001)}, who established a heterogeneous distribution in the mean values of air temperature for three stations evaluated in the Colombian Caribbean.

One of the relevant factors in the annual temperature cycle of the Colombian Caribbean is the Caribbean Low-Level Jet (CBNC) (^{Amador-Astúa, 1998}), which has been related to trade wind fluctuations, increased hurricane activity, and increased temperature in northern South America (^{Yepes et al., 2019}). This low-level jet in the Caribbean converges with the Northern Hemisphere summer (^{Amador-Astúa, 1998}; ^{Serna et al., 2018}) and with the temperature maxima reported in this study. In addition to this low-level jet, one of the main intra-annual modulators of climate in Colombia is the Intertropical Convergence Zone (ITCZ) (^{León et al., 2000}), which has a presence in the Colombian Caribbean, specifically during July. These two intra-annual scale phenomena coincide with the periods of one year and ~ 6 months present in the first and second harmonics calculated in our study from the cyclic descent analysis (Table II). This same analysis allowed us to identify periods between 3 and 7 years that are in agreement with ENSO oscillation periods (^{Rasmusson et al., 1990}), periods > 10 years that coincide with solar magnetic activity (^{Friis-Christensen and Lassen, 1991}, ¹⁹⁹⁴; ^{Connolly et al., 2021}; ^{Velasco-Herrera et al., 2021}). Both phenomena directly influence the climate on interannual scales, but there is evidence that during ENSO years, the CBNC and the ITCZ are affected (^{Builes-Jaramillo et al., 2023}). which reflects on air temperature in the Colombian Caribbean.

The evaluation of the effects of climate variability events on air temperature in the Colombian Caribbean was approached in this work by means of associations between time series of air temperature and the SOI, MEI, NAO, and Sunspot climate indices (Fig 3). The MEI was considered the most influential factor in air temperature fluctuations, with a statistical significance of 95% and a Spearman’s correlation coefficient average value of 0.37, followed by the SOI (r_sp = -0.31). The NAO did not present significant statistical associations, hence its effects are not considered relevant to air temperature fluctuations in the Colombian Caribbean. An entropy analysis between MEI and SOI performed by ^{Ortiz-Tánchez et al. (2002)} revealed different characteristics of both ENSO quantifiers, being the MEI index the one that shows quantitative evidence of mid-range correlations to the Southern Oscillation, since it presents a more constructive effect of its fluctuations in comparison with those of the SOI, which clearly established the onset of El Niño/La Niña and its effects beyond one month. In this sense, the MEI index behaves much better than the SOI in the description and predictability of ENSO, as we found in the study of its effects on air temperature in the Colombian Caribbean from the correlation coefficients and harmonics 3, 4 and 5. In addition to this, in the present study we were able to determine ENSO as the main influential factor in climate variability on an interannual scale, as has been reported in different regions of the world and some regions of Colombia (e.g., ^{Pabón-Caicedo et al., 2001}; ^{Rodríguez-Rubio, 2013}). Finally, although a significant influence of solar magnetic activity on climate variability has been reported (^{Buendía-Carrera et al., 2012}), for the Colombian Caribbean we found a low correlation statistically significant at 95%, indicating that although the influence of this phenomenon is not relevant, it is of vital importance to consider it in studies of climate teleconnections.

The evaluation of climate goes hand in hand with its tendencies, for which in the present study we determined for the Colombian Caribbean, based on non-parametric statistical techniques, that in the department of Cesar the average air temperature is increasing ~ 0. 01 ºC yr^-1 more than in the departments of Magdalena and Guajira (Fig. 5). ^{Carmona and Poveda (2014)} reported overall increases in mean air temperature in the Colombian Caribbean region after evaluating the trend of hydro-climatic variables throughout Colombia; they determined that 87% of the monthly minimum temperature series assessed presented increases between 0.01 and 0.08 ºC yr^-1. In this sense, according to the Third National Communication on Climate Change (TCNCC) of Colombia (^{IDEAM, 2017}), based on IPCC reports, the average temperature may increase ~ 1 ºC by 2040 with respect to the reference period 1976-2005 in the northern region of the country. The increasing trends found in this study may be related to the migration of inhabitants from rural regions to cities due to forced displacement, changes in land use, and desertification processes to which the Colombian Caribbean is being subjected. Therefore, in this region in particular, the growing trend could be related to an urban heat island effect, as mentioned by ^{Jones et al. (1990)}. In addition to overpopulation and changes in land use, increases in greenhouse gas emissions due to vehicle growth (^{Lamb et al., 2021}) could be another factor causing atmospheric warming and the increases in air temperature reported in this study.

Analyses of global trends conducted by trimesters, decades, or seasons indicate positive temperature trends associated with higher temperatures and therefore major impacts on the agricultural sector, which is considered by some as evidence of climate change (^{Al-Buhairi, 2010}; ^{El-Geziry, 2022}; ^{Mateus and Potito, 2022}; ^{De et al., 2023}). Finally, the presence of aerosols and particles of terrigenous origin from the Guajira deserts reported by ^{Orozco-Rodríguez (2021)} and from the Sahara Desert (^{Hernández, 2014}; ^{Poleo and Briceño, 2014}), could be affecting the energy balance in such a way that the reduction of incoming solar radiation to the surface (^{Slingo et al., 2006}) increases daily heating rates (^{Chen et al., 2017}), warming of the middle troposphere, and cooling of the lower troposphere (^{Mohalfi et al., 1998}), all of which could cause increases in air temperature in the Colombian Caribbean, as we report in the present study. Finally, the results obtained in this work on the spatial heterogeneity of mean air temperature values, the influence of the CBNC, the suitability of evaluating the effects of ENSO with the MEI index, and the positive trends, provide climate change scenarios that should be taken into account for the development of management plans aimed at thermal comfort, reduction of greenhouse gas emissions, and land use, thus favoring the environmental conditions of the Colombian Caribbean, its inhabitants, and the different species that coexist in it.

6. Conclusion

The analysis of the variability, cycles, and trends of mean air temperature in the Colombian Caribbean, based on data from IDEAM meteorological stations, allowed establishing that the reconstruction of time series by the NLPCA method guarantees the variability and mean values of air temperature in this region, such that the comparison metrics present mean values of 0.99, 1.8 × 10^-4 and 3.5 × 10^-7 ºC for Spearman’s correlation coefficient, mean absolute error, and bias. In addition, it was identified that the annual cycle of air temperature shows heterogeneity between the mean values and their spatial distribution, mainly due to the geographical position, orography, proximity to the sea, and humidity contributions of the continental environments belonging to this region.

The periodicity analysis established that the main modes of variability at the intra-annual scale are the CBNC and the ZCIT, while at a global scale ENSO predominated, which in turn showed better correlations than the other modes of climatic variability evaluated in order to identify the influence on air temperature fluctuations. This association analysis allowed us to characterize the MEI index as the most suitable for evaluating the effects of this mode of variability on air temperature.

Finally, the trends of the complete series oscillate between 0.01 and 0.07 ºC yr^-1, with a maximum in the department of Cesar. As for the monthly trends for the study period (1991-2020), it had a maximum between March and June in the department of Cesar (1.07 ºC yr^-1) and a minimum in La Guajira (0.2 ºC yr^-1) in March. These temperature increases are significant given the productive activities of crops and livestock in the region, which may be affected by them as well as by periods of drought and water resource availability. The results obtained in this study are very useful for decision-makers when considering the average air temperature variable, its variability, cycle, and trends, and therefore its effects on heat islands that currently exist in the population centers associated with the economic activities that take place in the region, thus allowing them to establish management plans aimed at thermal comfort and climate change mitigation measures.