Dynamics of meteorological time series on the base of ground measurements and retrospective data from MERRA‐2 for Poland

A comparison study has been performed to assess the dynamics of meteorological processes in Poland on the basis of meteorological time series of air pressure, air temperature and wind speed coming from 35 synoptic stations belonging to the Institute of Meteorology and Water Management—National Research Institute (IMGW‐PIB) and from the nearest grid points of the Modern‐Era Retrospective Analysis for Research and Applications, Version 2 (MERRA‐2) produced at NASA's Global Modelling and Assimilation Office from January 1, 2007 to October 31, 2016. Apart from comparative statistics, the differences in the multifractal properties of the time series were evaluated with the use of MultiFractal Detrended Fluctuation Analysis (MFDFA), both for hourly and daily data, showing a high degree of similarity between the MERRA‐2 and IMGW‐PIB series. For the air pressure and air temperature, not only were high determination coefficients (close to .99) between the time series coming from the two sources noticed, but there were also similarities with the MFDFA parameters. Lower correlations between the time series of the wind speed obtained from the two studied databases were observed, which was related to differences in the data structure and methodology of the measurements for specific IMGW‐PIB stations. Additionally, to verify data similarities coming from the IMGW‐PIB and MERRA‐2 databases, the correlations between specific multifractal parameters and the orography were estimated and compared. For the air pressure and temperature, a remarkably high correlation was found between the multifractal parameter α0 and the height above sea level of the measurement site. An analysis of the source of multifractality was performed, indicating that, for all studied meteorological elements and both data sources, the long‐range correlations prevail.


| INTRODUCTION
To observe the state of the atmosphere, an extended ground monitoring system of weather conditions has been created at a global scale of hundreds of thousands of meteorological stations, which collect enormous amounts of information each day concerning such fundamental weather parameters as: air temperature, air pressure, air humidity, precipitation, wind speed, and net radiation (Dahoui et al., 2020). The procedures for data collection have been standardized by the World Meteorological Organization (WMO, 2018). Long-term meteorological data (usually longer than 30 years) are often analysed for studies concerning climate change impacts on various aspects of human living conditions (Pirttioja et al., 2015;Fronzek et al., 2018;Ruiz-Ramos et al., 2018;Rodríguez et al., 2019). Existing historical observation data show spatio-temporal irregularities (Vose et al., 2012;Murat et al., 2016Murat et al., , 2018Stott et al., 2016;Gos et al., 2020) and were obtained with the use of a variety of instruments and measurement techniques which have been modified over time (Jones and Wigley, 2010). Additionally, inconsistencies in ground-level meteorological data are present in long time series due to issues with instrument relocations (Brohan et al., 2006), urbanization (Peterson et al., 1998), and land-use changes (Montandon et al., 2011). These ground-level data are often supplemented by measurements from satellites, balloons, aircraft, or ships, providing an opportunity for researchers to not only improve the spatial resolution of the data but to improve the global forecasting skill and quality of these longer time series given an improved understanding of the sources of inhomogeneity (Chahine et al., 2006).
An alternative to the direct measurements of the meteorological elements are retrospective analysis systems of weather conditions reaching far into the past, which use multiple sources of information (ground, satellite, and aerial data) and computer modelling (Simmons et al., 2007;Rienecker et al., 2011;Boilley and Wald, 2015). Reanalysis aims to assimilate historical atmospheric observational data spanning an extended period using a single consistent assimilation scheme to provide reprocessing of meteorological observations in a physically consistent manner. The output of the reanalysis is gridded datasets of studied meteorological elements for specific periods. Two main objectives for the creation of reanalysis weather datasets are to improve the available observational record for the early 20th century and to prepare datasets and assimilation tools for global reanalysis. The leading reanalysis systems are: MERRA and MERRA-2 from the National Aeronautics and Space Administration ; AgMERRA created within the Agricultural Model Intercomparison and Improvement Project ; CERA-20C; ERA-15, ERA-20C, ERA-40, and ERA-Interim from the European Centre for Medium-Range Weather Forecasts Dee et al., 2011;Laloyaux et al., 2018); NOAA-CIRES 20th Century Reanalysis V2c supported by the National Oceanic and Atmospheric Administration as well as the Cooperative Institute for Research in Environmental Sciences and the U.S. Department of Energy (Slivinski et al., 2019); North American Regional Reanalysis; Climate Forecast System Reanalysis; NCEP-DOE or NCEP-NCAR established at the National Centers for Environmental Prediction (Kalnay et al., 1996;Saha et al., 2010;Montealegre et al., 2014); and the Japanese 25-year Reanalysis (Kobayashi et al., 2015) from the Japan Meteorological Agency.
Reanalysis products have been successfully used for climate monitoring and atmospheric research, including the creation and verification of a high spatio-temporal resolution global precipitation dataset designed for hydrological modelling (Beck et al., 2017), capturing halfcentury trends in global air temperatures (Compo et al., 2013), evaluating the dynamics of the warming hiatus (Huang et al., 2017), modelling the response of plant diseases to climate change (Bebber et al., 2016), reconstruction of Atlantic tropical cyclone activity for the last millennium (Burn and Palmer, 2015), changes and interannual variability in daily scale temperature and precipitation extremes (Donat et al., 2016), analysis of the precipitation data properties for different land regions , selection of the regions vulnerable to the occurrence of extreme events (García-Marín et al., 2013;Baranowski et al., 2019), estimation of temperature-mortality associations (Royé et al., 2020), modelling of a global streamflow for the world's mediumto-large river basins (Candogan Yossef et al., 2017;Alfieri et al., 2020), analysing distributions of long-term pan evaporation across China , assessing cloud cover over the Arabian Peninsula (Yousef et al., 2020), assessing the impact of volcanic eruptions on climate extremes (Paik and Min, 2018), and the characterization of meteorological drought in India (Adarsh et al., 2019).
The reanalysis data quality must be evaluated as a fundamental issue before these data can be applied in multiple disciplines. There are many ways of accounting for uncertainty in the reanalysis data of meteorological elements, including comparison of data across various reanalysis sources and comparison with observational data in a specified time-scale (Compo et al., 2011;Alfieri et al., 2020). Such internal quantification of uncertainty enables many errors to be avoided in connection with inaccuracies in multisource data processing algorithms, including those for ground and satellite observations. Constant improvements and upgrades of reanalysis software also require verification of reliable local and regional observational data sets to improve their efficacy (Decker et al., 2012).
To assess the similarity of reanalysis data to observed values, a standard statistical approach is often chosen. The standard approach to detecting differences in the time series dynamics for meteorological elements is, for example, analysing the trends in the relevant meteorological values. This allows the confirmation of the spatial or temporal differences as a result of climatic conditions or changes within a specific period. This typical approach for time series is often not sensitive enough to predict differences in the time series coming from similar climatic zones or for series collected in the same regions using various types of instruments or methods. Additionally, this approach does not always consider that meteorological processes possess inherent features, such as being characterized by random fluctuations on different spatiotemporal scales and showing high nonlinearity. The dynamics of meteorological elements can be evaluated using multifractal analysis (Kantelhardt et al., 2002), which enables the detection of long-range correlations in noisy, nonstationary time series of respective variables . The multifractality of meteorological elements has already been extensively studied, indicating strongly self-similar properties in most climate variables, including precipitation (de Lima and de Lima, 2009;Gemmer et al., 2010;Baranowski et al., 2015), ozone concentration (Jiménez-Hornero et al., 2009;Pavon-Dominguez et al., 2013), air and surface temperature (Yuan et al., 2012;Jiang et al., 2013;Burgueño et al., 2014), or wind speed (Krzyszczak et al., 2017a;Laib et al., 2018). Additionally, the spatio-temporal variabilities of multifractal parameters of meteorological elements have been recognized at various scales (Hoffmann et al., 2017;Krzyszczak et al., 2017b). An interesting but insufficiently investigated issue in this topic is preserving self-similarity in reanalysis data as compared to the multifractal properties of those ground-level data measurements .
The aim of our study was to compare the dynamics of meteorological time series derived from MERRA-2 (the second Modern-Era Retrospective analysis for Research and Applications) and ground-based data for Poland obtained from the Institute of Meteorology and Water Management-National Research Institute (IMGW-PIB) database. The specific aims were focused on comparisons of: • basic statistics, including the Mean Error (ME), Root Mean Square Error (RMSE), and coefficient of determination (R 2 ) of the MERRA-2 and ground series data; • multifractal properties of MERRA-2 and ground-based decadal meteorological time series for different locations in Poland; and • dominant sources of multifractality in the meteorological time series through the analysis of original, shuffled, and surrogated time series. A temperate climate ('Dfb' according to the Köppen-Geiger classification) characterizes Poland, with relatively warm summers (where the average temperature ranges between 18 and 30 C depending on the region) and cold winters (with average temperatures around 3 C in the northwest and −6 C in the northeast; Wibig and Glowicki, 2002). Poland is strongly influenced by oceanic air from the west, cold polar air from Scandinavia and Russia, as well as warmer, sub-tropical air from the south. The weather of the area is highly unpredictable and varied. The northern coast, overlooking the Baltic Sea, and the western part of the country are influenced by the milder oceanic climate (Cfb).

| Data from synoptic weather stations
The analysed ground-based time series came from the IMGW-PIB database in Poland (https://danepubliczne. imgw.pl); its observation and measurement system covers the entire territory of Poland and consists of 62 hydrological and meteorological stations. The IMGW implements ratified conventions with the World Meteorological Organization, the International Civil Aviation Organization, and the Safety of Life at Sea, and it cooperates with European and global meteorological and hydrologicalmeteorological services.
The ground-based datasets contained hourly time series of three meteorological quantities: air pressure, air temperature (at 2-m height), and wind speed (usually at 10-m height, but it varies from 9 to 20 m). The analysed data covered the period from January 1, 2007 to October 31, 2016 and came from 35 meteorological stations that evenly covered the whole country. The criterion for choosing the 35 stations from among the 62 synoptic stations in Poland was that, within the entire 10-year period, they contained no disturbed hourly data. Some of the synoptic stations that were not included conducted measurements only three times a day or contained long gaps (more than 1 month), and the localisation of some of the excluded stations changed during the studied period. Daily data from the three studied meteorological variables, obtained by averaging the hourly data, were also analysed in the study.

| MERRA-2 data
An overview of the MERRA-2 system can be found in Gelaro et al. (2017). MERRA-2 was generated with the use of the Goddard Earth Observing System-5, Data Assimilation System, and Gridpoint Statistical Interpolation that uses an incremental update procedure every 6 hr (Wu et al., 2002;Kleist et al., 2009). The system was supplied with numerous additions compared to the features of the previous MERRA system, including newer microwave sounders and hyperspectral infrared radiance instruments, as well as additional data types (McCarty et al., 2016).
The analysed MERRA-2 grid data (284 points) were obtained from NASA (http://disc.sci.gsfc.nasa.gov/mdisc/) servers with a resolution of 0.5 north-south and 0.625 east-west. The MERRA-2 dataset contained hourly and daily time series of three meteorological elements: air pressure and temperature at 2-m height and wind speed at 10-m above the ground for each grid point. The analysed data covered the period from January 1, 2007 to October 31, 2016, the same as for the IMGW-PIB data. The studied datasets consisted of continuous data sequences averaged over the indicated interval and time stamped with the central time of the interval (Bosilovich et al., 2015).
The map of all of the locations of the ground stations (35 points drawn in green) and the MERRA-2 nearest land-based grid points (indicated with red crosses) used in the analyses are presented in Figure 1. Sea-surface grid points were excluded from the study, as our preliminary analyses indicated that meteorological elements at these points possess different statistical properties than those that are land-based.
All the maps in this article are represented by the geographical longitude and latitude expressed in degrees. The description of the reference ground points of the IMGW-PIB stations together with the coordinates of the nearest MERRA-2 grid points are presented in Table 1.

| Basic comparative statistics of ground-based and MERRA-2 time series
The first step in comparing the raw meteorological data from two sources, IMGW-PIB and MERRA-2, was to calculate the basic statistics, including the mean error (ME), the mean squared error (RMSE), and the coefficient of determination (R 2 ). ME and RMSE were defined according to the following formulas: F I G U R E 1 Comparison of MERRA-2 grid points (red crosses) and the localization of ground stations (green dots) used for calculation T A B L E 1 Locations of the meteorological ground stations and the respective MERRA-2 grid data points used for the comparison of the where N is the number of measurements, and e i = X i − Y i is the difference between the measured value X i of a given meteorological element in a specific IMGW-PIB point at a specific time and the value of Y i , representing the respective measurement from the MERRA-2 database. The ME is the average value of e i . The RMSE is the square root of the average squared values of the differences between the measurements from IMGW-PIB and the respective measurements from MERRA-2 for a given meteorological element. These errors have the same units as their respective meteorological elements (air pressure: kPa, air temperature: C, and wind speed: mÁs −1 ). The coefficient of determination R 2 was used to explain how much of the variability could be caused by its relationship to the respective values of the element from the MERRA-2 datasets. R 2 was calculated by the formula: To calculate these statistics, the RStudio integrated development environment for R version 0.97.551 was used (R Core Team, 2018).
Additionally, to compare the results of the IMGW-PIB and MERRA-2 observations for the three considered meteorological elements, a scatter plot with the colour indicating the data density was created in MATLAB software version R2018b (Sanchez, 2020).

| Comparison of the multifractal properties of the IMGW-PIB and MERRA-2 time series
To compare the dynamics of the time series from two different sources (IMGW-PIB and MERRA-2), their multifractality was evaluated with the use of MultiFractal Detrended Fluctuation Analysis (MFDFA), which was elaborated on by Kantelhardt et al. (2002) and is described in detail in many other articles (Kalamaras et al., 2019;Krzyszczak et al., 2019;Liu et al., 2019). The MFDFA for nonstationary time series is performed using five steps. Firstly, the noises are converted into random walks by subtracting the mean value from each actual value within the time series and by integrating the time series to determine the associated profile. Then, each profile is divided into non-overlapping segments of equal lengths, s. This procedure is repeated starting from the opposite ends of the series. Another step is fitting a proper polynomial function by a least-squares for each of the created segments and subsequently determining the variance. All segments are averaged to obtain the qthorder fluctuation function F q (s): where N s is the number of non-overlapping segments, and ν = 1…2N s is the number of segments. The scaling behaviour of the functions is determined by analysing log-log plots F q (s) versus s for each value of q. For multifractal time series, F q (s) increases (for large values of s) as a power law F q (s) $ sh(q) with the generalized Hurst exponent h(q) depending on q. The multifractal spectrum is obtained using the relationship τ(q) = qh(q) − 1 and using the Legendre transform f (α) = qα-τ(q), where α = dτ/dq.
The schematic representation of a multifractal spectrum is shown in Figure 2. In this article, the following parameters were analysed: • α 0 : the value of α at the maximum of the power spectrum f(α); • a s : the asymmetry parameter, which can result in negative or positive values for a left-or right-skewed shape, respectively; • w: the width of the spectrum (the difference between α max , the fractal dimension representing the mildest processes, and α min , the fractal dimension corresponding to the most intense processes in the analysed series).
When the values of the α 0 parameter become lower, the underlying process loses its fine structure and becomes more regular in appearance. The width parameter, w, is responsible for indicating the richness of the signal structure; when it increases, the length of the range of fractal exponents in the signal becomes higher, suggesting more developed multifractality. The symmetry parameter a s provides information about the structure of the signal and the possible occurrence of extreme events (in the case of strongly negative values of this parameter). A right-skewed spectrum indicates strongly weighted high fractal exponents and provides information about fine structures of the signal (Telesca and Lovallo, 2011).
To determine which source of multifractality in the original time series is dominant, transformations of longrange correlations or a broad probability density function were performed. These transformations included shuffling by random permutations (RP), Amplitude Adjusted Fourier Transform (AAFT), and Fast Fourier Transformation (FFT), as suggested by (Theiler et al., 1992). According to Kantelhardt et al. (2002), the random shuffling of the elements of the time series removes the long-range correlations, but the multifractality due to the broadness of the probability density function cannot be removed with a shuffling procedure. In general, if surrogate data obtained either by AAFT or FFT lead to multifractals not differing from the respective multifractals of the original data, the broad probability density function of the time series does not influence the multifractality. Additionally, an FFT without the amplitude adjustment procedure destroys only nonlinear correlations in the analysed signal, leaving the linear ones. Therefore, using this transformation, it is possible to distinguish whether the nonlinear or linear long-range correlations are dominant within long-range correlations. If the nonlinear correlations of variability prevail, the FT spectrum of a multifractal time series is very narrow (Nagarajan, 2007;Hall, 2014;Mali, 2014;Baranowski et al., 2015).

| RESULTS AND DISCUSSION
4.1 | Results of the statistical comparison between ground-based data and MERRA-2 data There is a good correlation between the IMGW-PIB and MERRA-2 data, especially for air pressure and air temperature. A much lower correlation exists for both the hourly and daily wind speed data. This is caused by the precision of the wind speed data delivered to the users (IMGW-PIB delivers data without decimals, whereas MERRA-2 delivers data with two decimals). Additionally, this issue may be connected to the difference in the altitude of the wind speed measurements between IMGW-PIB (9-20 m above the surface level for some meteorological stations) and MERRA-2 (10 m above the surface level). For the other two meteorological elements, there is no difference in the height of the measurements and the accuracy of the data presentation. In Tables 2 and 3, the basic statistics are presented for the three studied elements on the basis of evaluating the correlation between the IMGW-PIB and MERRA-2 data coming from 35 particular locations. This means that, for each location, the R 2 , RMSE, and ME were calculated for individual elements, F I G U R E 2 Schematic representation of the main parameters of a multifractal spectrum D 0 , information carrier; D 1 , information dimension; α, singularity strength; q, order of the fluctuation, a s , asymmetry parameter T A B L E 2 Comparative statistics of hourly data from IMGW-PIB and MERRA-2 for three studied meteorological elements

Statistical parameters
Air pressure IMGW-PIB and MERRA-2 Air temperature IMGW-PIB and MERRA-2 Wind speed IMGW-PIB and MERRA-2 Another way of comparing the IMGW-PIB and MERRA-2 data was done via presentation on the density scatter plot (Figure 3a,c,e for hourly data and Figure 3b, d,f for daily data). The R 2 , RMSE, and ME values presented in this figure differ from those obtained in Tables 2  and 3 because they refer to all the pairs of values within a 10-year period at 35 locations (3,017,282 measurements). Figure 3 also presents the distribution of the measuring points' density expressed as normalized values for the studied meteorological elements. The scatter density plots in Figure 3 indicate that the highest density of the points for air pressure occurred in the range of 97-102 kPa, that for air temperature occurred in a range from −5 to 20 C, and that for wind speed occurred in the range of 0 to 4 mÁs −1 . These are the most frequently occurring values of these meteorological elements in the climatic conditions of Poland. In Figure 3, linear regression lines (green colour) and 1:1 lines between the IMGW-PIB and MERRA-2 data (red colour) are presented.
Accordingly, Tables 2 and 3 and Figure 3 indicate that, for air pressure, there is nearly no difference in the R 2 , RMSE, and ME values between hourly and daily data. For air temperature, the daily data have slightly higher R 2 values, lower RMSE errors, and the same ME errors. The lowest correlations between the IMGW-PIB and MERRA-2 data and the highest differences between the hourly and daily data were observed for the wind speed. In Figure 3e (hourly data), a low R 2 = .56 and a relatively high RMSE = 1.94 mÁs −1 indicates a very low correlation, which slightly improves for the daily data (R 2 = .70 and RMSE = 1.56 mÁs −1 ). For the daily wind speed data, the regression line is above the identity line for the daily data, and the ME value is negative, which means that MERRA-2 gives slightly higher values of wind speed than IMGW-PIB reports. In the case of hourly data, the ME value is also negative, but the regression line indicates that MERRA-2 overestimates the wind speed for low values and underestimates it for high values. The analysis for individual quarters shows that no differentiation of synoptic contexts can be identified for air pressure and wind speed when comparing the IMGW-PIB and MERRA-2 data. On average, MERRA-2 slightly underestimates the air pressure values, but this underestimation was on the same level for all the seasons for both hourly and daily data. In contrast, the wind speed values were overestimated, with a similar deviation from the values measured for individual quarters, regardless of whether hourly or daily data were used. The R 2 and RMSE were on the same level for each quarter for the air pressure and wind speed data. Different behaviour and influence of synoptic context was identified for the air temperature data. Under a regime of strong synoptic control (winter), we observed that the data from MERRA-2 were underestimated, whereas, under a regime of weaker synoptic control (summer), they were overestimated. The correlation and precision of the data fit was on the same level for each quarter (R 2 and RMSE had similar values for T A B L E 3 Comparative statistics of daily data from IMGW-PIB and MERRA-2 for three studied meteorological elements

| Comparison of multifractal properties of IMGW-PIB and MERRA meteorological time series using MFDFA
In this study, we decided to perform MFDFA analysis to confirm and extend the results of the basic comparative statistics of the time series coming from IMGW-PIB and MERRA-2 presented in the previous section. Since the purpose of the study was not only to compare meteorological data obtained from two separate sources but also to compare hourly and daily data, multiple segment sizes, s, were used in the analysis, which would reflect similar time lags being considered in MFDFA for these two aggregations. After some experimenting, it was decided that, for hourly data s min = 120 and s max = 9,000, whereas, for daily data, s min = 5 and s max = 375. This F I G U R E 3 Scatter density plot of three studied meteorological elements derived from IMGW-PIB and MERRA-2 databases. Hourly data are presented in a,c,e and daily data in b,d,f means that the minimal segment size for both hourly and daily data covered 5 days, and the maximum segment size covered 375 days. Choosing the minimum and maximum segment size is a particularly important part of MFDFA analysis because the qth-order fluctuation function F q (s) depends on it. During selection of the multiple segment sizes, Kantelhardt (2009) recommended diminishing discrepancies or an abrupt decline in the scaling exponents. Additionally, an analysis similar to the one presented in Palus (2008) was conducted to show the behaviour of the scaling exponent τ(q) for different values of q. Palus (2008) used synthetic multifractals F I G U R E 4 Comparison of similarity spectra of hourly and daily meteorological time series which did not influence the long-range correlations, broad tails in the probability density function, nor the unprocessed time series. By analysing a broad range of qϵ (−10;10), they indicated that scaling exponents for FT and AAFT surrogates are similar in the range of qϵ(−3;6), otherwise they should differ considerably depending on the source of multifractality. In our study, the range of q was limited to [−5;5] to diminish the distortions connected with the so-called 'freezing' effect (Kantelhardt et al., 2006). According to the methodology of Palus (2008), we analysed surrogate (AAFT and FFT), shuffled (RP), and original data. We also experimented with selecting the optimal m value by running MFDFA multiple times with different m values (ranging from 1 to 4). Finally, we determined that the behaviour of F q (s) versus s was most stable for m = 1. High similarity exists between respective parameters of multifractal spectra coming from IMGW-PIB and MERRA-2 time series both for hourly and daily data (Figure 4). The left side of Figure 4 presents a comparison of hourly data, and the right side presents the daily data. In each plot, the red and green colours refer to the IMGW-PIB and MERRA-2 data, respectively. The lines show the medians, and the shaded areas indicate the 5th-95th percentiles of the regions, covering similarity spectra derived from the 35 studied locations. For air pressure and air temperature, the shaded areas and median lines of IMGW-PIB and MERRA-2 overlap, indicating high degrees of similarity in the underlying processes.
With respect to the wind speed, the right sides of the areas covering the spectra for IMGW-PIB are broader and are moved to the right in relation to the areas of the MERRA-2 spectra. The time aggregation process strongly influenced the parameters of the spectra. The daily spectra were usually broader than the respective hourly spectra for the three studied meteorological elements (higher width of the spectra). The daily time series spectra of the studied elements had higher asymmetry than did the hourly ones. The least significant differences existed between the values of the α 0 parameter between the spectra of the hourly and daily time series. These results are confirmed by the parameters of the multifractal spectra of the original time series, presented in Table 4 for hourly data and Table 5 for daily data. In each cell of these tables, the upper number is the mean value of a parameter (from the average of 35 spectra at multiple locations), whereas the lower numbers are the minimum and maximum values of this parameter. The respective multifractal parameters of the studied meteorological elements coming from the IMGW-PIB and MERRA-2 databases indicate strong similarities (mean values and range of change). It is also evident from these tables that the highest degree of similarity between the hourly and daily multifractal parameters occurs for α 0 (for all studied meteorological elements).
In the case of air pressure and air temperature, α 0 was significantly higher than for the wind speed, both for hourly and daily data. A lower α 0 for wind speed means that time series of this variable are characterized by better correlated underlying processes, loss of subtle structure, and, thus, greater regularity over the course of a given time series. For hourly and daily data, the asymmetry parameter a s of the air pressure is always positive, which means that the spectrum is right-skewed, denoting strongly weighted high fractal exponents, which correspond to the subtle structure of the series (time series are characterized by a multifractal structure, which is insensitive to local fluctuations with large amplitudes). Similar results were found for the daily air temperature spectra, whereas, for other elements, this parameter takes on both negative and positive values. Negative values of a s (leftskewed spectrum) for some stations suggest the existence of complex structures at the level of large fluctuation amplitudes in the respective series; they also suggest that the time series are characterized by a multifractal structure, which is insensitive to local fluctuations with small amplitudes. This means that a low weight can be attributed to low fractal exponents, which in turn suggests the frequent occurrence of extreme events (Telesca and Lovallo, 2011).
It is interesting that, for the air temperature, the hourly a s data take on negative and positive values, while the daily data take on only positive values. That result suggests that some information about the structure of the series is lost when aggregating the hourly data to daily F I G U R E 5 Relation between α 0 parameter and the altitude for hourly (left side) and daily (right side) data of air pressure and air temperature data, which is in agreement with the results obtained by Krzyszczak et al. (2017b). We also noticed that the width of the spectrum w was lowest in the case of hourly wind speed data, which means that the respective time series of this variable had smaller variety in the individual fractals and lower complexity in the signal structure (less 'rich' multifractality) than the other analysed variables did. A more detailed analysis of the differences in multifractal properties of time series from MERRA-2 for Polish territory is presented in Baranowski et al. (2019).
Because the time series of the studied meteorological elements came from the entire territory of Poland, we decided to check whether the altitude of the stations has any impact on the multifractal properties of the analysed time series. Three parameters of the multifractal spectra were included in this analysis; however, the width w of the spectra and the asymmetry a s did not have high correlations with the altitude. Only the air pressure and the air temperature α 0 values indicated relatively high correlations with the altitude (with their absolute values ranging from 0.63 to 0.82), both for hourly and daily data, which are presented in Figure 5. In this figure, various altitude ranges are indicated with separate colours. The linear regression equations and their coefficients of determination (R 2 ) are presented in the plots separately for IMGW-PIB and MERRA-2 data. It results from Figure 5 shows that the α 0 parameter increases for the air pressure and decreases for the air temperature with an increase in altitude. When comparing the relationship between α 0 and altitude for the data coming from IMGW-PIB and MERRA-2, it was observed that-both for the hourly and daily time series-the fitted linear regression lines have nearly the same equations for air pressure, but, for the air temperature from the MERRA-2 data, the α 0 values move upward nearly parallel to the IMGW-PIB α 0 data. This confirms that the altitudinal effect on the structure of the air temperature time series obtained from these two data sources is different.
It has been suggested in previous articles (Burgueño et al., 2014;Laib et al., 2018;Baranowski et al., 2019) that the terrain orography has an impact on the multifractal spectra of meteorological elements (daily and hourly air temperatures, wind speed, wind direction, and air pressure). The complex topography of the Swiss territory influenced the differentiation of multifractal properties of daily wind speed (Laib et al., 2018). In Poland, there is no such high differentiation in elevation as there is in Switzerland. Poland covers mainly lowlands and highlands, with altitudes not exceeding 300-400 m, except for in a small area in the southern region with mountains with  average peaks of $1,300-1,400 m and not higher than 2,500 m a.s.l. In our study, we had too few measuring points with altitudes exceeding 300 m a.s.l. to confirm the Laib et al. (2018) results. Additionally, to acknowledge the source of multifractality for the studied time series of the three meteorological elements, spectra of the shuffled and surrogated time series were created. It was found by Kantelhardt et al. (2002) that the artificially transformed time series was able to distinguish between the dominance of the multifractality due to a broad probability density function for the values of the time series and the multifractality due to different long-range correlations for the small and large fluctuations. The randomly shuffled series (obtained by generating random permutations of the array elements of the time series) make it possible to remove any temporal correlations. If the spectra narrow significantly (width parameter decreases) after shuffling, long-term correlations are predominant in the multifractality of the data. The multifractality coming from broad distributions can be tested using an AAFT (Theiler et al., 1992). If the multifractality in the time series is due to a broad probability density function, the spectra obtained for the AAFT surrogate data are much narrower than for the original ones (Min et al., 2013;Mali, 2014). Additionally, for nonlinearity tests of the time series, FFT was used. Fourier transform surrogates without the amplitude adjustment procedure destroy nonlinear correlations in the analysed signal, leaving only the linear correlations. Therefore, it enables us to distinguish whether the nonlinear or linear long-range correlations are dominant within the overall long-range correlations. If the nonlinear correlations of variability prevail, the FFT spectrum of a multifractal time series is very narrow (similar to a monofractal spectrum). If the linear correlations of variability prevail, the FFT spectrum is very rich (high width), and we can say that the process is strongly nonlinear. The results of the analysis of the sources of multifractality within the studied time series of meteorological elements are presented in Table 6 for hourly data and in Table 7 for daily data. These tables are composed of three parts: spectra from AAFT surrogates, spectra from FFT surrogates, and RP spectra (shuffled data).
The parameters of multifractal spectra of surrogated data (AAFT and FFT) are very similar to the respective parameters of the decomposed original data for both hourly and daily data. Accordingly, parameters of the shuffled data differ significantly from the parameters of the spectra obtained from the original and surrogated data. The shuffled spectra of all the studied meteorological elements are characterized by the α 0 parameter being T A B L E 7 Mean, min and max values of MFDFA spectra parameters for studied meteorological elements (daily data)  (Livina et al., 2003(Livina et al., , 2011 and wind loads (Kareem and McCullough, 2013). In Figures 6 and 7, the averaged data (from 35 localisations) of the scaling exponents τ(q) for the three studied meteorological elements are presented as functions of q. In these figures, the FT, AAFT, and RP surrogates are compared with the original data. The solid lines show the mean values, and the bars show the standard deviations of 100 randomly generated realizations of these transforms.
For high values of q, both negative and positive, the multifractality of FFT and AAFT differ from the multifractality of the original data. This effect is visible F I G U R E 7 The averaged scaling function τ(q) taken from 35 sites for original daily decomposed data (black line) and the mean of 100 realizations of their AAFT (green solid line), FFT (red solid line) and RP (blue solid line) surrogate data. Vertical bars indicate 95% of distribution of appropriate surrogates for all analysed meteorological elements, both for the IMGW-PIB and MERRA-2 plots. On the other hand, a quite different pattern of the scaling exponent change with q was observed for the RP series, compared to the original and surrogate series. Observed straight lines for all the studied meteorological elements of hourly data suggest that the multifractality was destroyed by shuffling, and typical monofractals occurred. This is clearly seen for the hourly data, but, for the daily data, the RP curves are slightly bent for all three meteorological elements. For both the hourly and daily wind speed, the RP curve more closely matches the AAFT and FFT curves, which suggests that the source of multifractality can be more complex in the case of this meteorological element. The change in the scaling exponents with q for the wind speed differs from the results obtained by Baranowski et al. (2015), who noticed larger differences among the AAFT, FFT, and original data. This discrepancy is probably related to the different times spans for these two studies (30 years and 10 years) as well as the averaging of 35 meteorological data points from various locations.

| CONCLUSIONS
A comparison between meteorological time series derived from MERRA-2 (the second Modern-Era Retrospective analysis for Research and Applications) and ground based data for Poland obtained from IMGW-PIB database was made using statistical evaluation and multifractal analysis. Those methods are complementary to each other, as they inform us about different properties of the time series. The standard approach based on statistical analysis determines the averaged state and the size of fluctuations around this state or estimates the trend and seasonality of the time series of a given agrometeorological variables (Balling et al., 1998). In turn, the MF-DFA method allows to assess inherent structure (system state at remote points that affects the local evolution of the system) of these time series.
High similarity between the MERRA-2 and IMGW-PIB multifractal properties of the meteorological time series was found, especially for the air pressure and temperature data. This similarity was confirmed not only by high determination coefficients ($.99) between time series coming from these two databases but also the resemblance of the averaged parameters from the MFDFA. The ranges of the change in the multifractal parameters (α 0 , a s , w) are comparable for the IMGW-PIB and MERRA-2 datasets for both hourly and daily data. Lower correlations between the time series of the wind speed obtained from the two studied databases probably come from the differences in the measurement heights and differences in data structure (IMGW-PIB delivers wind speed data without decimals).
For the air pressure and air temperature, a high correlation was found between the multifractal parameter α 0 and the height above sea level of the measurement site. This result suggests that, for air pressure and air temperature, the dynamics of the processes change significantly with altitude, which should be precisely determined in cases of interchangeable use of the data from various sources.
Undoubtedly, many practical implications of the similarities between the IMGW-PIB and MERRA-2 time series datasets can be derived. Firstly, any existing gaps in ground-based meteorological time series from specific locations can be filled successfully by MERRA-2 data, which possess a similar structure and properly reflect the dynamics of the original series. This is relevant for both present-day processes as well as those from 30 years ago due to the reanalysis made within the MERRA-2 database. The resolution of the MERRA-2 data grid is also satisfactory to deliver valuable data for supplementing information from any of the ground stations. The distances between the ground stations and MERRA-2 grid points differed from 3.9 to 62.2 km in this study, but the similarities between the two sources of meteorological information seemed to be satisfactory for many applications, such as the modelling of crop growth and development under various climatic conditions, as suggested previously by Baranowski et al. (2019). The results suggest that MERRA-2 data can be especially useful for analyses of climate dynamics changes and can be used interchangeably with the groundbased data.