1. Introduction

The chaos theory assumes the disordered, disorganized, confused, incoherent, in dynamic systems (^{Edwards and Penney, 2009}; ^{Ramírez, 2010}; ^{Fernández, 2016}; ^{Lestayo, 2021}). The Andean atmospheric system is chaotic, limiting knowledge about its behavior (^{Richter et al., 2005}; ^{Feo et al., 2009}), which becomes a problem given its influence on daily human activities, hence the need to find alternatives for its description.

It is necessary to determine the variability of atmospheric dynamics on a small scale in the Andean region, specifically in the city of Riobamba, given the interest in its particular conditions such as height above sea level and equatorial geographic position, when concluding air and sea currents that influence the variability of climatic behavior on a small to micro scale.

The structure of the Earth’s atmosphere is extremely complex, made up of layers and sublayers that vary at different scales. However, due to technological development, it has not been possible to describe it in its entirety, several theories have emerged, such as the chaos theory, which perfectly adjusts to it and allows us to study certain particular characteristics (^{Motter and Campbell, 2013}; ^{Pino et al., 2018}). It is in the troposphere where atmospheric turbulence, caused by thermal or mechanical imbalances, produces an increase in the chaotic movements of the air, thus favoring the phenomenon of diffusion and transport of energy and matter, in which all human activities take place and therefore it is necessary to know its behavior (^{Ivancevic and Ivancevic, 2007}; ^{Haro et al., 2016}; ^{Pino et al., 2018}).

The principles of conservation of mass, momentum, and energy are fundamental in describing both the internal processes of the atmosphere and their external interactions. While the gases that constitute it endure significant variations in temperatures and pressures, about the height above sea level (^{Almendra et al., 2022}). These dynamic tropospheric processes, caused by mechanical turbulence related to wind and friction with the solid surface, as well as convective turbulence produced by the flow of heat between the ground and the air, result in heating of the ground during daylight hours (^{Rezaei et al., 2022})

For the present investigation, the procedures formulated by van Ulden and Hostlang associated with the study area and the surface energy thermal fluxes are adopted. These determine the variation of superficial atmospheric dynamics or atmospheric stability according to the variation of the air’s physical properties. Riobamba, located between 1º38’3”-1º4’ S and 78º39’- 78º40’36” W in a plateau at 2750 masl, surrounded by high mountains (^{Haro et al., 2016}).

The dynamic processes of the atmosphere that surround it result in turbulent systems with low descriptive and predictive precision using nonlinear equations of the dynamics that cannot be described analytically. However, there are models developed in mid-latitude countries that fulfill this function with a certain error bias and do not adjust to the local conditions in the city of Riobamba, given the dependence on semi-empirical parameters to be adjusted (^{Haro et al., 2016}). This is the reason for studying the small-scale dynamic variation of micrometeorological parameters, appealing to known models of atmospheric dynamics and adjusting them to local conditions. Thus, the surface heat flux, the latent heat flux, the friction velocity and the Obukhov length were determined, knowing the behavior of atmospheric microdynamics.

The study recalls the chaos theory, which studies complex systems, sensitive to initial conditions, where small variations imply considerable differences in future behavior and make long-term prediction impossible. Such a phenomenon occurs (^{Bolotin et al., 2017}; ^{Fowler and McGuinnes, 2019}) in systems of deterministic rigor, that is, its behavior can be completely determined by knowing its initial conditions and computerizing approaches and equations of chaos theory in given atmospheric problems (^{Medvinsky et al., 2017}; ^{Weixuan et al., 2022}).

The chaos theory is an effective tool to study complex or chaotic systems, since in real systems it is common to find signals that appear to be non-random behavior. Characterized by high sensitivity to initial conditions and unpredictability over time, and through variables in space where the dimensions represent dynamic variables (^{Hegger et al., 1999}; ^{Setoudeh et al., 2022}).

Few theories are adapted to the system’s characteristics we are dealing with, and one of the best adapted to complex atmospheric systems is the chaos theory, which has demonstrated its efficiency in various systems with similar characteristics, and in recent years it has gained great importance thanks to the development of computers. Throughout time, the atmosphere has dominated life on the earth’s surface and all species have had to adapt to its conditions. Despite scientific and technological development, man is no exception, and knowing his behavior allows him to make better use of his resources and prevent natural disasters due to atmospheric phenomena.

It is important to apply these principles that have been developed in other latitudes to different systems for validating the results, previous adjustments, and conditions of the area. Following the theoretical and practical procedures of the theory, it has been possible to achieve specific results for the study area which contribute to the knowledge of the atmosphere that adjusts to similar systems and allows us to quantify the variations that occur in it by means of parameters.

2. Materials and methods

The present investigation declares a longitudinal design given the time of sampling in four years; quantitative approach, by taking and analyzing data from weather stations, mediating mathematical models. According to the level of depth in the object of study, of a descriptive and explanatory type, declaring scientific methods such as the analysis and synthesis of informational sources, mathematical statistical method, inductive-deductive method.

Regarding the collection of data from meteorological stations (Fig. 1), variables such as temperature, wind speed, humidity, solar radiation, and atmospheric pressure are assumed in hourly averages. The date was recorded at the weather stations using sensors, and any missing data was imputed using linear interpolation.

For data processing, the calculation of surface energy fluxes in the atmospheric boundary layer is based on ^{Beljaars and Holtslag (1999)}, and de Bruin and Holtslag (1982). Those that are produced by the interactions between the Earth’s surface and the atmosphere, those that generate the convection of air in the atmosphere, making it highly turbulent and unpredictable under certain conditions. In addition, the method according to ^{van Ulden and Holtslag (1985)} is adopted, associated with the study of surface flows in the atmospheric boundary layer by changing the properties and physical characteristics of the circulating air mass on different surfaces.

The basic parameters that describe air conditions, such as temperature and humidity, are transformed from the contact surface and propagated to higher layers. Defined as a boundary layer, where energy is transferred from the surface to the atmosphere and vice versa (de ^{Bruin and Holtslag, 1982}) through turbulent flows of sensible (Q _h ), latent (Q _g ) and momentum heat. In this way, ^{Kiely (1999)} indicates that the total energy balance of the volume element based on the surface that contains vegetation cover and the surrounding atmosphere will be used to produce evapotranspiration. It must be evaluated, translating the result into units of evaporated water according to the simplified energy balance equation:

The terms R _N and Q _g are measured with the pyrometer using the relationship of ^{Kasten and Czepelak (1980)} to determine the cloudiness with Q _g as a function of R _N (Eq. [2]; de ^{Bruin and Holtslag [1982]}) with small thermocouples buried in the ground. Sometimes the heat flux stored in the soil (C _a ) is also neglected, and this causes important errors. Therefore, neglecting this term will produce more error, the shorter the measurement interval, where on occasions C _a is estimated empirically as a fixed function of R _N (^{US-EPA, 2004}). Then Q _h and Q _L cannot be measured separately, but through the quotient between both flows (Q _h /L _g ), called the Bowen relation. Starting from a surface energy balance, it is expressed in Eqs. (2) and (3):

where a = 0.1 for rural areas and a = 0.3 for urban areas in this study (^{Kasten and Czepelak, 1980}), β is a constant equal to 20 W m-2. and α = 0.75 for urban environments. The parameter S is defined by Eq. (4):

Considering the kinetic energy variation equation for turbulence, static fluctuations are established from accessible variables, a turbulent energy equation for closed surfaces is represented in Eq. (5)

where the term on the left represents the total energy flow, while the first term on the right represents the thermal fluctuation of energy, the second the mechanical fluctuation of energy, the third the dissipation of energy due to the effect of friction, and the fourth the energy flow term that combines the turbulent transport term and the pressure term (^{Sang et al., 2009}). Now it is possible to determine the Obukhov length (L), which is a certain height in which the thermal production is balanced with the mechanics (Eq. [6):

The atmospheric similarity theory states that ∂U/∂z = u _* /kz, where is the friction velocity and k is the von Karman constant. By applying the Monin-Obukhov similarity theory, the friction speed u _* is determined as a function of the wind speed U _z and the height z (^{Teixeira et al., 2008}; ^{Sang et al., 2009}) (Eq. [7]):

With the von Karman constant k (k = 0.4), z ₀ length of surface roughness of the place that is determined by tables (^{Marrero, 2011}).

Substituting Eq. (6) and considering Z = L, the length of Obukhov length (Eq. [8]) is calculated as:

where the parameter on the right is considered surface and can be determined indirectly through meteorological measurements (wind speed, temperature, cloudiness, etc.). Equation that is generally used as the atmosphere’s dynamic state indicator (^{McNaughton, 2009}; ^{Rodríguez et al., 2015}).

The stability function ψ _M is determined by L < 0 (unstable), as indicated by Eq. (9):

From Eqs. (8) and (9), the friction velocity and the Obukhov length can be found by interaction, initially taking ψ _M = 0 and L = ∞ until reaching an approximate variation of 1% of two successive values of L for values greater and less than zero and with the same criteria, for the stable case L > 0 the relationship is applied (de Bruni and Holtslag, 1982; ^{Doll et al., 1985}):

where

The data characterized from the chaos theory begins with the reconstruction of the series in the space of the phases and thus determine the delay time and the fit dimension. These allow to represent said data in a multidimensional space and to study its dynamic properties such as the fractal dimension and the entropy of the system. The important coefficients, such as Lyapunov, determine if the system is chaotic or not, in addition to predicting the future behavior of the data (^{Motter and Campbell, 2013}; ^{Abderrahim et al., 2014}).

Regarding the time series, results in a sequence of observations {S _n = S(x _n )}, which is generally scaled {S _n }, and by itself, does not adequately represent the (multidimensional) phase space of dynamical systems. Then, the technique to determine the multidimensional structure using the available data, or the delay method, must be employed. As vectors exist in a new space of the lace space, values of the lagged time series are formed, obtaining the following matrix:

Where S _n is the reconstructed matrix of the data series, τ the delay time, and m the dimension of fit.

Thus, the reasonable delay time is determined, a dimensionless integer value that allows the series to be reconstructed as a matrix to be represented in a multidimensional space, recalling the mutual information method from ^{Fraser and Swinney (1986)}. Unlike the autocorrelation function, the mutual information also considers the non-linear function, calculated with Eq. (12) (^{Haro et al., 2016}).

where S _m is the mutual information function, p _i is the probability of finding a time series value in the i-th interval, and p _ij is the joint probability that an observation is found in the i-th interval and later in the j-th interval. The expression calculation is executed if the delay value τ finds a minimum, it is a good option to use it as delay time. In general, the first minimum of Eq. (12) is assumed as an acceptable value of the delay time (^{González, 2004}; ^{Pino et al., 2014}).

The number m of elements or fit dimension (Takens’ theorem) (^{Sancho, 2016}), establishes that if {S _n } is the sequence of measurements of a dynamical system, the fitting dimension with proper lag time provides a unique picture of the original ensemble if m is large enough. This measurement is applied to scalar measurements N and sets of lace vectors m with dimension N - τ(m - 1).

These must be considered for the calculation of average quantities in phase space. From literature on the choice of optimal parameters m the method for determining the minimum size m is proposed, where the dimension is accepted when the function is zero (^{González-Miranda, 2004}; ^{Haro et al., 2016}). The idea is that for each point in the time series, the nearest neighbors must be found in m-dimensional space. Calculating the distance s→i-s→j between the points and calculating the ratio given in Eq. (13):

If R _i exceeds a threshold value R _t this point is called a false neighbor, the criterion is that the dimension is accepted when the fraction of the points R _i > R _t is small enough or zero. Chaos arises from the exponential growth of infinitesimal disturbances. To guarantee the analysis of this instability, there are Lyapunov exponents, which measure how much two trajectories move away (^{Bochi and Rams, 2016}; ^{Barreira, 2017}):

where r _i (t) is the semi-major axis of an ellipse, r the small initial radius, and t the large time, in general we get:

λ ₁ ≥ λ ₂ ··· ≥ λ _d

Variables that allow characterizing an attractor, as follows:

From this result, the Kolmogorov-Sinai entropy can be defined in Eq. (15).

Or the dimension of Kaplan-Yorke (1979) according to ^{Stolz et al. (2017}):

with j-positive Lyapunov coefficients.

3. Results

The data is processed from weather stations over four years at one-hour intervals using software designed by the authors of this research These stations are located in the canton Riobamba, which is in the province of Chimborazo, Ecuador. Specifically, they are situated in the Ecuadorian highlands at an altitude of 2700 masl. The three stations used in this study are equipped with Vaisala sensors, and their characteristics are presented in Table I.

Basic micrometeorological parameters such as the Obukhov length, surface heat flux, and latent heat flux are calculated based on van Ulden’s proposal determining atmospheric dynamics. To calculate the variation of these parameters, the freeware TISIAN is used. Determining the delay time, fit dimension, Lyapunov coefficients, Kaplan-Yorke dimension, Kolmogorov-Sinai entropy, and thus obtaining the following results:

Taking the data from the three stations, initially, noise reduction was carried out because certain processed parameters lacked data correlation, particularly the friction velocity, which ultimately is not presented. However, the other parameters, which are presented in Figure 2, could be processed.

The delay times values for each of the variables are not very different between weather stations or parameters, as can be seen in Figure 2. Figure 3 illustrates the strange attractors of surface heat in the three stations. It can be observed that there are not many differences between them.

The superficial heat flux produced by the air-soil interaction generates superficial thermal fluxes, which rise rapidly during daylight hours, decreasing at night, with a high variability during the year and its attractors have irregularities as shown in the figure.

A higher variability is observed regarding latent heat, especially between QUIMIAG and the other two stations (Fig. 4).

The latent heat flux, which plays a crucial role in determining atmospheric thermal conditions, exhibits significant diurnal variations. It tends to decrease during nighttime and shows a variable behavior throughout the day. The reconstructed attractors in two dimensions display irregularities (Fig. 5), but the data from QUIMIAG station exhibit slight differences.

The Obukhov longitude (Fig. 6) fluctuates significantly due to sudden changes in temperature, particularly during the early and late hours of the day. The attractor exhibits a sufficient level of regularity among the three stations, except for minor deviations observed at the ESPOCH station, where a few values exceed 1000 m of fluctuation.

4. Conclusions