Time-Series Analysis of Isotope Composition of Precipitation in Zagreb, Croatia

Abstract

Air temperature and precipitation data (1976–2021), stable isotope composition (δ¹⁸O, δ²H) data, and deuterium excess (1980–2021) data were analyzed using principal component analysis (PCA), Fourier analysis (FA), and wavelet analysis (WA). The PCA represented each month by a single dot in the diagram, and month 1 and month 7 were clearly distinguished. The FA and WA gave the 12-month period for all parameters, but the strongest power was for temperature, then δ¹⁸O and δ²H, and finally for the precipitation amount and deuterium excess. Both Pearson’s r and Spearman’s ρ correlation coefficients gave similar values for δ²H—δ¹⁸O and temperature—δ²H, δ¹⁸O correlations.

1. Introduction

A statistical technique that deals with time-series data is called a time-series analysis. Time-series data mean that the data are in a series of particular time periods or any other intervals. It is used by many industries in order to extract meaningful statistics, characteristics, and insights. Among time-series analysis, Fourier analysis (FA) and wavelet analysis (WA) have proven to be powerful tools for detecting temporal variations in time-series data.

Principal component analysis (PCA) is a dimensionality-reducing method that is used to reduce the dimensionality of large data sets by transforming a large set of variables into a smaller one that still contains most of the information in the large set. Such a small dimensionality set can be easily visualized and analyzed.

The strength of the relationship between the two variables is called a correlation coefficient, and its values range between −1.0 and 1.0. A correlation of −1.0 is a perfect negative correlation, while a correlation of 1.0 is a perfect positive correlation. A correlation of 0.0 shows no linear relationship between the two variables.

Few stations have enough data to perform wavelet data analysis. About one-sixth of the number of complete cycles can be confidently found [1]. Therefore, at least 6 years of the data record is needed for accurate and precise determination of 1-year seasonal fluctuations. The longer the observation time, the higher the accuracy and precision in evidencing the cycles and their time span.

Most of the WA were used to study global long-term periods [2,3,4,5,6,7]. Some were used to teleconnect large-scale climate indices and hydrochemical and isotopic characteristics of a karst spring [8] or to study spatio-temporal variability of rainfall and streamflows over the last decades [9].

Only a few studies have been applied to tritium analyses [10,11] or stable isotopes in precipitation [12]. Wavelet analysis was used to determine tritium variation in precipitation for several stations with long-term data and correlate them with solar activity fluctuations [10] and determine solar cycle periodicity in tritium activity in precipitation in Zagreb from 1996 to 2019 [11]. Time-series analysis (both FA and WA) was applied to a 7-year record (2012–2018) of δ¹⁸O, δ²H, deuterium excess in precipitation, as well as temperature and precipitation amount [12].

The history of monitoring isotope data in precipitation at Zagreb is much longer than the minimum of seven years—tritium has been monitored since 1976, while stable isotope data are available for 1980–2006 and from 2012 onward. The series of data is rather scarce in Europe [13], and it deserves to be analyzed in many different ways. Statistical analysis of ³H activity concentration, δ²H, δ¹⁸O, and deuterium excess and comparison with basic meteorological data (temperature, amount of precipitation) was presented in [14], but no principal component analysis was performed, and no Fourier and wavelet analysis was performed. Mean annual temperature showed a statistically significant increase of 0.07 °C per year. Annual mean δ¹⁸O and δ²H values showed an increase of 0.017‰ and 0.14‰ per year, respectively, resembling an increase in the temperature. The mean yearly precipitation amount did not show any statistically significant change, but the annual values showed higher dispersion from the mean during the last 20 years. The distribution of the monthly amount of precipitation moved to the second half of the year, with the maximum in September. The mean deuterium excess remained constant over the years, but a change of distribution was observed: a decrease in deuterium excess in the first half of the year and an increase in the second half (maximum in November) due to air masses coming from the eastern Mediterranean.

Long-term (1980–2021) isotope data (δ²H, δ¹⁸O, deuterium excess) of monthly precipitation in Zagreb were studied here using statistical time-series analyses. Data for meteorological parameters for 1976–2021 (mean monthly air temperature and the monthly amount of precipitation) were also included. We used the data presented in [14] complemented with the new data from 2019–2021 (Supplementary File, Table S1). Principal component analysis, Fourier analysis, and wavelet analysis were used to determine what kind of results could be obtained. The principal question was whether this kind of analysis could show any dependence on temperature. Finally, we compared the two correlation coefficients (Pearson’s r and Spearman’s ρ) to determine any possible difference between them.

2. Materials and Methods

Two frequency analysis methods were used: periodogram (Fourier analysis) and wavelet. The algorithm for both methods used supports evaluating unevenly spaced data points.

Unevenly and missing spaced data points could be converted to evenly spaced points in a few different ways, such as interpolation, setting missing data points to zeros, and others. However, practice showed that these methods are not reassuring. Therefore, for periodogram analysis, we chose the Lomb periodogram [15]. In our case, N data points $(h_i\equiv h\left(t_i\right), i=1,\dots ,N)$ are measured at unevenly sampled times $t_i$. The Lomb normalized periodogram is given by

$$P_N\left(\omega \right)=\frac{1}{2{\sigma }^2}\left\{\frac{{\left[{{\displaystyle \sum }}_j\left(h_j-\overline{h}\right)\cos \omega \left(t_j-\tau \right)\right]}^2}{{{\displaystyle \sum }}_j{\cos }^2\omega {\left(t_j-\tau \right)}^2}+\frac{{\left[{{\displaystyle \sum }}_j\left(h_j-\overline{h}\right)\sin \omega \left(t_j-\tau \right)\right]}^2}{{{\displaystyle \sum }}_j{\sin }^2\omega \left(t_j-\tau \right)}\right\}$$

(1)

where $\tau$, $\overline{h},$ and ${\sigma }^2$ are defined by

$$\tan \left(2\omega \tau o\right)=\frac{{{\displaystyle \sum }}_j\sin 2\omega t_j}{\cos 2\omega t_j}$$

(2)

$$\overline{h}=\frac{1}{N}{\displaystyle \sum }_1^N{h}_i$$

(3)

$${\sigma }^2=\frac{1}{N-1}{\displaystyle \sum }_1^N\left(h_i-{\overline{h}}^2\right)$$

(4)

and angular frequency $\omega =2\pi f,$ f being the frequency. The significance level of any peak in $P_N\left(\omega \right)$ is given by $1-{\left(1-e^{-P_N\left({\omega }_0\right)}\right)}^M$ for some concrete $\omega$₀, and $M$ stands for the number of independent frequencies. The method used for fast calculation of $P_N\left(\omega \right)$ and $M$ can be found in [16].

To analyze signal transient and/or non-sinusoidal nature, wavelet analysis was chosen. The algorithm for unevenly spaced data proposed by [17] was used. For the mother wavelet, Morlet wavelet $f\left(z\right)=e^{-c{z}^2}\left(e^{iz}-e^{-1/4c}\right)$ with c = 1/72 was used [1].

Our data can be viewed as a multidimensional point consisting of five variable vectors $\left(x_1,x_2,x_3,x_4,x_5\right)$, where each point corresponds to a particular month of the year. The vector components $x_i$ corresponds to temperature, δ¹⁸O, δ²H, precipitation, and deuterium excess. Such arranged data are suitable for principal component analysis (PCA) [18] in the points where all variables are present—in our case, 416 data points in total. Two principal components (PC1 and PC2) were used to present the data point, and the confidence ellipse and centroid were calculated. The confidence ellipses were calculated assuming bivariate Gaussian distribution. The centroid point $\left(X,Y\right)$ for $k$ points $\left(x_i, y_i\right), i=1,\dots ,k$, was derived as follows:

$$X=\frac{1}{k}{{\displaystyle \sum }}_{i=1}^k{x}_i$$

(5)

$$Y=\frac{1}{k}{{\displaystyle \sum }}_{i=1}^k{y}_i$$

(6)

Both Pearson’s and Spearman’s methods were used to measure the correlation between the data. Pearson’s correlation coefficient r evaluates the linear relationship between two continuous variables, while the Spearman correlation coefficient ρ is based on the ranked values for each variable rather than on the raw data. Both the correlation coefficients give a positive (or negative) correlation value, but the values can be different because Pearson’s r measures δ, the linear relationship between the variables, while Spearman’s ρ only measures monotonic relationships.

New stable isotope data (δ²H and δ¹⁸O of Zagreb’s precipitation), related to the period 2019–2021, were determined at the Laboratory for Spectroscopy of the Faculty of Mining, Geology, and Petroleum Engineering, University of Zagreb, with a Liquid Water Isotope Analyzer (LWIA-45-EP, Los Gatos Research, San Jose, CA, USA). Data were analyzed by the Laboratory Information Management System (LIMS) [19]. The measurement precision of duplicates was ±0.19‰ for δ¹⁸O and ±0.9‰ for δ²H.

3. Results

The temperature in 2019–2021 showed an increase in the lowest temperature and steady highest temperature in the last 5 years [13], which gave 2019 the title of the hottest year of the whole period of 1976–2021. Similar conclusions are also valid for δ¹⁸O and δ²H data. The annual mean deuterium excess is constant throughout 1980–2018, but with higher values in autumn than in spring, and this difference increased over time [14]. Deuterium excess for 2019–2021 was around 8‰ for January–August, while for September–December, it was higher than in all other periods, including 2012–2018 [14], reaching the average value of 14‰ in November. The average amount of precipitation did not change in 1976–2018, but larger fluctuations around the mean values were observed during the most recent period [14]. Precipitation in 2019–2021 was average (± 1 standard deviation), equal to the average value of 1976–2018.

3.1. PCA Analysis

PCA analysis (Figure 1) shows data points for months 1 and 7, together with their confidence ellipses. Those two months, each representing a different season, can be statistically separated in the PC1-PC2 space. The PC1 axis accounts for 56.27% of the variability, and PC2 for 25.04% of the variability in the data. We can also see that temperature, δ¹⁸O, and δ²H are more closely aligned with the PC1 axis, and precipitation and deuterium excess are linked with the PC2 axis. This leads to the conclusion that the winter/summer difference is most pronounced in temperature, δ¹⁸O, and δ²H changes.

Additionally, Figure 1 shows centroids from all twelve months, with their path in the PC1-PC2 space denoted here by vectors marked with the numbers of the month. Such a representation further emphasizes that warmer months are mostly on the right side of the PC1 axis, corresponding to a higher temperature, δ¹⁸O, and δ²H. On the PC2 axis, precipitation and deuterium excess dominate. Months 3, 4, and 5 are separated from 11, 10, and 9, suggesting that the latter have more precipitation and a greater excess of deuterium. Here, the difference is less pronounced than on the PC1 axis.

3.2. FA Analysis

Fourier analysis (FA) of the temperature, δ¹⁸O, and δ²H values (Figure 2a,b and Figure S1) showed the strongest peaks. All peaks were well above the significance level of p = 0.05. The peak appeared at 12.01 months since those values exhibit strong annual seasonality. However, the increase in temperature due to global climate changes in recent decades cannot be seen [14].

FA of the precipitation amount (Figure 2c) and deuterium excess (Figure 2d) showed more frequencies, although the mean frequency is still one year, and others are not statistically significant. The 1-year peak is still well above the limit of p = 0.05, but the power of the peak (~30) is significantly lower than that of temperature (~250) or δ¹⁸O and δ²H (~120).

3.3. WA Analysis

The temperature shows a typical graph of seasonal fluctuation within a period of 12 months (Figure 3a), localized at the (12 ± 1) month period. δ¹⁸O (Figure 3b) and δ²H (Figure S2) graphs also show typical seasonal graphs with slightly lower localization in a 12-month period. Here, we should mention that δ²H and δ¹⁸O graphs were obtained for 1980–2006 and separately for 2012–2021 because of a large gap in data from 2007 to 2011. However, the same dominant period is observed in both cases.

The WA graph for precipitation and deuterium excess (Figure 3c,d) showed an average period of 12 months, but with much larger fluctuations compared to δ²H and δ¹⁸O—they contained more non-dominant periods. The deuterium excess graph is also divided into 1980–2006 and 2012–2021 periods with the same characteristics.

WA analysis confirms that the dominant 12-month periodicity is stationary in time.

3.4. Correlations

Correlations among the studied quantities were quantified by both Pearson’s (r) and Spearman’s (ρ) correlation coefficients for linear and monotonic correlation, respectively (

Table 1. Pearson’s r and Spearman’s ρ correlation coefficients.

ρ r	δ18O (‰)	δ2H (‰)	Precipitation (mm)	Temperature (°C)	Deuterium Excess (‰)
δ18O (‰)	-	0.986	0.084	0.799	−0.238
δ2H (‰)	0.988	-	0.117	0.799	−0.098
Precipitation (mm)	0.089	0.127	-	0.246	0.225
Temperature (oC)	0.793	0.799	0.235	-	−0.137
Deuterium excess (‰)	−0.277	−0.126	0.235	−0.155	-

Numbers in bold—significant correlation.

and Figure 4). Both gave similar results, and the strongest correlations were found between δ²H and δ¹⁸O (both r and ρ values > 0.98), temperature and δ²H (both r and ρ values > 0.79), and temperature and δ¹⁸O (>0.79), justifying the linear and monotonic correlations. The other quantities showed much weaker or no correlation at all.

4. Conclusions

The time series analyses justified previously observed behavior of temperature, precipitation amount, and stable isotope data in precipitation concerning seasonal periodicity. All used statistical methods indicate a similar behavior pattern of temperature, δ¹⁸O, and δ²H. Both FA and WA gave the highest power periodicity for the temperature, followed by δ¹⁸O and δ²H, and finally, the precipitation amount and deuterium excess. FA showed a peak at 12 months, which was confirmed by WA, suggesting a stationary periodicity in time.

However, this kind of analysis does not give information on the long-term temporal change of the parameters, such as the dependence on temperature. This could be an important parameter for relating the stable isotope data in precipitation with the current climate change.

Thus, the results of the current paper present complementary results to the previously published ones. In the future, any analysis of precipitation data for at least 7 years of isotope data will not be considered complete without the following: principal component analysis, Fourier or wavelet analyses, δ¹⁸O vs. δ²H relation, and δ¹⁸O vs. temperature or precipitation amount.