During few last years, climate change including global warming which is attributed to human activities and also its long-term adverse effects on the planet's functions have been identified as the most challenging discussion topics which have arisen many concerns and efforts to find the possible solutions. Since the warmth arising from Earth's landscapes affects the world's weather and climate patterns, we decided to study the changes in the Land Surface Temperature (LST) patterns in different seasons through non-linear methods. Here, we particularly want to estimate the non-integer dimension and fractal structure of the land surface temperature. For this study, the (LST) data has been obtained during the daytime by the (MODIS) on NASA's Terra satellite. Depending on what time of the year data has been collected, temperatures change in different ranges. Since equatorial regions remain warm, and Antarctica and Greenland remain cold, and also because altitude affects temperature, we selected Riley County in the U.S. state of Kansas, which does not belong to any of this type locations and we are interested to observe the seasonal changes in temperature in this county. The results of the present study show that the Land Surface Temperature (LST) belongs to the class of fractal process with non-integer dimension. With increasing temperature around the world, photosynthesis and microbial and plant respiration and as a result releasing carbon dioxide to the atmosphere also increase [26]. During recent decades, climate changes and its long-term impacts such as accelerating ice melt in poles, warming oceans, and rising sea level and many other aspects of life entire earth has been widely recognized by many different organizations and researchers [5,18]. According to the most recent Intergovernmental Panel on Climate Change (IPCC) report (IPCC, 2014 [38]) greenhouse gases (GHG) are at the highest level which caused the most significant changes in earth temperature in history [18]. Human activities have significantly raised the concentration of greenhouse gases in atmosphere and caused global warming which is affecting adversely all life of the earth [38,44]. Climate change and global warming result in human activities by changing the Land Surface Temperature (LST) make significant changes in ecosystem structure and consequently function [43,50]. For example, they cause acceleration in the environmental status changes on the Tibetan Plateau (Zhong et al. ( 2011)) [50]. Many 2/21 studies have also documented the effect of land-use and land-cover change which are caused by human activities on land-atmosphere interactions (Boisier et al., 2013, Boisier et al., 2012, de Noblet-Ducoudr e et al., 2012) [12,13,17]. States at Risk which is a project to demonstrated the effects of climate changes in all 50 states of United States, divides the climate changes impacts into five categories: extreme heat, drought, wildfires, coastal flooding and inland flooding. There are many real-world examples of these five impacts of climate changes in different parts of the earth. For example, the devastating brush-fires happened in Australia in 2019 because of high temperature is one the most dramatics recent disasters result in climate changes which followed by extreme heat and long drought beside killing people and animals (Source(s): SafeHome.org, States at Risk). Across the United States, we can see different impacts of climate changes in daily life such as extreme heat during summer and intolerable cold during winter, destructive storms, wildfires prolonged bush-fires in plains and many more that have a significant impact on economy and other aspects of lives. (LST) plays a key role in environmental studies and climate researches [18]. The (LST) demonstrates the temperature of any surface of the Earth in a specific location. Since, a MODIS satellite collects the data through atmosphere to the ground, this surface can be snow, ice, grass, leaves of trees, roof of buildings. Therefore (LST) is not the same as the air temperature that we hear in the daily weather report. (LST) is a combination of vegetation and bare soil temperatures and because these two factors vary very fast in response to any solar radiation and or aerosol variations, (LST) demonstrates rapid changes as well. Scientists by collecting this data want to find the relationship between increasing atmospheric greenhouse gases and land surface temperature, and generally speaking, discovering the impact of increasing (LST) on glaciers, ice sheets, permafrost, and the vegetation in Earth's ecosystems (Source(s): (MODIS)). In the United States, southwest including central plains in Kansas confront unprecedented drought. Increasing the earth temperatures and reduction in rainfall will cause more extreme drought in the future than before effects adversely the urban dwellers and also farmers across this state. Using the data -winter 2022. Climate is a complex system since not all the variables in this system can be observed and because of uncertainty, our knowledge to predict is limited. To possibly describe complex natural processes, Benoit Mandelbrot in 1983 introduced fractals as immensely useful object in numerous scientific fields [34]. Mathematically speaking, fractals are infinitely irregular self similar objects which are characterized by power law or scaling exponent over a wide range of scales [33,34]. Thanks to Mandelbrot' research in fractal geometry, earth scientists are able to make probability and predict the size, location, and time of future natural disasters by measuring past events. If a certain process can be characterized by a single scaling exponent, then it is called mono-fractal. However, if a phenomena requires a large number of scaling exponents (can be infinite) to characterize their scaling structure, then it is called multi-fractal [27,31,42]. In studying real world time series data, we may deal with databases which exhibits power law or self similarity a wide range of size and time scales. For this type phenomena, we need to apply non-linear techniques to characterize the complexity and self similarity of data because the traditional time series techniques fail to classify databases with a variety range of scaling features [5-8, 30, 35]. The present study investigates the self-similarity and power-law structures of the Land Surface Temperature (LST) time series data using some non-linear techniques and fractal geometry.
We start with statistical regression analysis on the Land Surface Temperature (LST) databases of a selected region, the Riley County in the U.S. state of Kansas. This database has been taken using Google Earth Engine, the (MODIS) on NASA's Terra satellite. Because of the irregularity in the nature of the (LST) database in different days of the year, we use the non-linear methods such as power-law analysis, fractal dimension analysis and multifractal spectrum technique to characterize the structure of data. The vibration analysis using power spectral densities (PSD) method has been performed to find whether some type of power-law scaling exists for various statistical moments at different scales of these databases. Then, we use Discrete Wavelet Transform (DWT) and Wavelet Leader Multifractal (WLM) analysis to discover the possibility that these time series databases belong to the class of multifractal process for which a large number of scaling exponents are required to characterize their scaling structures. Finally, a non-linear analysis called the Fractal Dimension (FD) analysis using Higuchi algorithm has been carried out to quantify the fractal complexity of data. The study area is Riley County which is located in the U.S. state of Kansas. According to the U.S. Census Bureau, its total area is 622 square miles (1610km Data: (LST) of Riley County for February, April, July, and October (2021) The present study investigates the non-linear patterns of daytime (LST) databases which have been collected using Google Earth Engine and over Riley County in the U.S. state of Kansas. Google Earth Engine is a combination of multi-petabyte catalog of satellite imagery and Geo-spatial databases with planetary-scale analysis capabilities [20]. Google Earth Engine helps to track changes and or measuring the differences on the Earth's surface and has been used widely to find map trends. To explore and quantify the trends and patterns in (LST) database, we perform a time series regression analysis. Time series regression is a statistical technique which uses experimental data to predict the future (LST) based on the current (LST) responses by transferring the dynamics of our predictors which is the height. This forecasting method provides a great understanding of the structure of the dynamic system which is (LST) database [21]. We start with defining the design matrix Ht of current and past height data and ordered bt time and we call this matrix as regressor matrix. Next, we use ordinary least squares (OLS) method for the linear regression (LR) equation: where  is the linear estimator parameter, Lt represent the (LST) database. The model (1) helps to find the linearity between the (LST) variable Lt and regressor matrix Ht. If we replace the linear part of equation ( 1) with a non-linear function f (Ht, t), we can better understand the relationship between the (LST) variable Lt and regressor matrix Ht as the following form: See figure (2). When we are working with stationary databases with constant frequencies throughout the time intervals of data, Fourier transform helps to visualize the data over time, however, this technique fails when we have datasets which are changing in time in response to a particular stimulus, i.e. they have time varying features [15,22]. The reason is the main limitation of Fourier transform which gives a broad power spectrum by integrating the frequency of time series data along the whole data. Because our Days (LST) databases are non-stationary, we will apply time-frequency methods such as Continuous Wavelet Transform (CWT). Wavelet-based approaches provide the ability to account for temporal (or spatial) variability in spectral character. Although wavelet analysis is relatively new compared with Fourier analysis, its use has become widespread in recent years, and so the theory will be described only briefly here [16]. Mallat (1998) or Strang and Nguyen (1996) contain excellent and detailed summaries of wavelet analysis. The continuous wavelet transform (CWT) of a dataset h(t) is given by (Mallat, 1998) [16,32] CWT(u, s) (3) where s is scale, u is displacement,  is the mother wavelet used, and * means complex conjugate. The CWT is therefore a convolution of the data with scaled version of the mother wavelet. Of course, the time coordinate t in equation ( 3) could equally well be the spatial coordinate x if profile data were being analyzed. Similarity measures are becoming increasingly commonly used in comparison of multiple datasets from various sources. Semblance filtering compares two datasets on the basis of their phase, as a function of frequency. Semblance analysis based on the Fourier transform ET(f) [11] where h(t) is dataset, f represents the frequency and t is the time, suffers from problems associated with that transform, in particular its assumption that the frequency content of the data must not change with time (for time-series data) or location (for data measured as a function of position). To overcome these problems, semblance is calculated here using the continuous wavelet transform. When calculated in this way, semblance analysis allows the local phase relationships between the two datasets to be studied as a function of both scale (or wavelength) and time [16,45]. When we define the Fourier transform (4) of time series datasets h1(t) and h2(t), the semblance S or the difference in the phase angle at each frequency can be obtained using the following formula s = cos(f) [14,47] where R1(f ) and R2(f This is a complex quantity and its amplitude which is called the cross-wavelet power defined as with local phase  arctan Im(CWT a,b ) Re(CWT a,b ) The local phase   (, ) measures the phase correlation between two time series data. Now, we define the wavelet-based semblance as the form Here n > 0 is odd and integer. For S = 1, our time series datasets are inversely correlated, for S = 0 there are un-correlated, and for S = +1, they are correlated. However, the wavelet-based semblance (8) fails to demonstrate the correlation between two noisy time series data because of of the lack of amplitude information [16].
Therefore, we define the vector dot product of the two complex wavelet vectors at each point in scale and position [16] D = cos n  |CWT a x CWT b * | (9) 9/21 As we can see from equation ( 9), D is a combination of the wavelet-based semblance (8) and the cross-wavelet power (7), which helps to find the phase correlations of the larger amplitude components of the noisy dataset. In figure (7), we have plotted the Dot product D of Land Surface Temperature (LST) of Riley County (figure (4)) during February, April, July and October (2021) using the formula (9) where dark red corresponds to large-amplitude signals with semblance equals +1, and yellow corresponds to large-amplitude signals with positive semblance between 0 and +1. The semblance results show that (LST) dataset and height datasets are highly correlated. Power spectral density (PSD) analysis is one of the methods in frequency-domain analysis which has been used frequently to analyze the signal vibration. Power defines the magnitude of the PSD which is the mean square value of the given time series data. PSD as a function of frequency demonstrates the distribution of a time series data over a spectrum of frequencies and its magnitude is normalized to a single Hertz bandwidth. PSD curves display the continuous probability density function of each measure, and it is obtained via multiplying each frequency bin of fast Fourier transform (FFT) by its complex conjugate to get a real spectrum. Because, the nature of (LST) database is non-linear and non-stationary, the welch (PSD) method with overlapped segmentation which is an averaging estimator technique has been applied to study the complex fluctuations in (LST) structures. Then, we use least square approximation technique to fit a linear regression model to the logarithm of power spectral density of (LST) database. Finally, we compute the slope for each regression model of (LST) data. In figure (8) In fractal processes, we can find a scaling relationship between power and frequency f in the spectral domain. The PSD results in figure (8) represent the fractal structure in (LST) database via a linear, negative slope of fitted least square lines. As result, the (LST) time series data cannot be obtained by one or a finite set of subsystems, and different components in each fractal process act at different time scales. Even though, PSD revealed the power law structure in (LST) databases, however, it failed to classify these four groups of data. Multi-fractal Analysis and Discrete Wavelet Transform (DWT) We have shown that complex structures in (LST) databases in figure (4) exhibit long-range correlation, i.e. they are self-similar: statistically, vibrations at one time scale are similar to those at multiple other time scales. To discover whether these four groups of (LST) data belong to the class of multi-fractal processes or they are mono-fractal, we have plotted the scaling exponents of the Land Surface Temperature (LST) of Riley County during February, April, July and October (2021) in figure (9). The non linear relationships between scaling exponents and the statistical moments could be a sign of multi-fractality in the structure of data, however, we need to check this claim via multi-fractal analysis. Multi-fractal analysis helps us to understand weather the (LST) databases are multi-fractal, i.e. means that a larger number of scaling exponents are required to characterize the scaling structures of time series data. To investigate scaling law of a multi-fractal process, the standard partition multi-fractal function was first introduced which was based on generalized fluctuation analysis. However, it failed to fully characterize the scaling behavior of the non-stationary time series data [10]. Arneodo et al., 1995b [4], Bacry et al., 1993 [9], Muzy et al., 1993Muzy et al., , 1994 [36] [36], developed a statistical techniques based on the Continuous Wavelet Transform (CWT) to characterize the multi-fractal description of singular measures. Arneodo [49]). In this study, we focus on multi-fractal analysis using (DWT) and (WL) [23,24,28,40,41,48,49]. The wavelet analysis associates the dimension of the fractal sets to Holder exponent H( ) to quantify the spectrum of singularity of the pointwise regular function F [28]. The Holder exponent of a fractal process F ( ) can be defined as follows: Definition 0.1. [41] A fractal process F ( ) satisfies a Holder condition, when there exist H( ) > 0, such that We can find H( ) for constant F from the coarse Holder exponents as |1 2 |< The following sets may be defined to extract the geometry of a signal with varying d, these sets describe the local regularity of signal. We call the map which is a compact form of the singularity structure of the fractal process F , the multi-fractal spectrum of F [41]. In a global setting, to describe the complexity of a signal, we may need to count the intervals over which the fractal process F evolves with Ho lder exponent H( ) and it gives an estimation of dim(E [d] ). The grain exponent which is a discrete approximation to h( ) can be written as the following form [41] 12/21 k Therefore, the grain multi-fractal spectrum has the form [23,24,40] log where From multi-fractal analysis in figure (10), we can easily see a clear loss of multi-fractality in all four groups of (LST) time series data, i.e. they are homogeneous and mono-fractal since their spectrum displays a narrow width of scaling exponent. One of the key concepts in theory of fractal geometry is the concept of fractal dimension. There are a variety of techniques to approximate the dimension of an irregular object. For example, the dimension of an object can be obtained via covering it with small boxes. The Hausdorff-Besicovitch dimension or box-counting dimension is one of the most often used methods to estimate the fractal dimension which works in this way [39]. To find box-counting dimension, we consider a square L x L and partition it into grids of linear length .
If N = (L/) 2 defines the total number of boxes, then, the total measure would be [19]: where we denote the measuring unit by  =   and try out a number of different . If we choose  = 1 then M   as   0. This means that the length of a square is infinite, which makes sense. If we try  = 3, then M  0 as   0. This means that the volume of a square is 0, which is also correct. For  = 2, we can find the true dimension of square. Using this argument, we define the Hausdorff-Besicovitch dimension as follows: first, we need -covering measure, which is the summation of the total measure 13/21 m m m m in all the "boxes" denoted by Vi, i = 1, 2, . . . , subject to the condition that the union of the boxes covers the object E, while the size of the largest Vi is not greater than .. Although, this method can successfully measure the self-similarity of a fractal process, however, it fails when the sudden changes happen in the irregular time series data sets [37]. To solve this problem, a variety of different non-linear methods such as Higuchi algorithm, power spectrum analysis, and Katz algorithm have been developed by different researchers [25,29]. To find the fractal dimension of each (LST) timeseries data, we apply the Higuchi Algorithm [25]. At first, we define a finite time series Finally, we use the least-squares approximation to estimate the slope of each regressed line as Higuchi fractal dimension for an optimal value of time interval k = 500 (optimal means there is no change in Higuchi fractal dimension after this value). We have estimated the fractal dimension of the Land Surface Temperature (LST) of Riley County during February, April, July and October (2021) and plotted their regression models for each time series in figure (11). From the results in figure (11), we can compare the fractal dimension of different (LST) data as follows FD February < FD July < FD October < FD April Although, fractal dimension using the Higuchi algorithm is a good index for comparing the self-similarity and complexity in fractal structure of different (LST) data sets, however, it fails to classify four groups of (LST) data sets. During recent decades, the rapid modification in land's surface by human activities has caused numerous environmental impacts. Land surface plays an important role in the structure of climate system and many related organizations put so much efforts to build observational land surface systems. Since controlling the Land Surface Temperature (LST), strongly impact the climate change in local and global scale, we decided to analyze the complexity of (LST) database to better understanding its pattern and 14/21 To explore the differences in the fractal patterns and complex dynamics of different (LST) data sets, we established our analysis based on different quantitative and analytical non-linear methods to find a computational framework in classification of (LST) databases. We started with vibration analysis using power spectral densities (PSD) method to estimate the exponent from realizations of these processes and to find out if the data of interest exhibits a power-law PSD. We performed multifractal analysis using discrete wavelet methods to discover whether some type of power-law scaling exists for various statistical moments at different scales of (LST) databases and also to explore if our database follows multifractal structure for which a large number of scaling exponents are required to characterize their scaling structures. Finally, we computed the Higuchi fractal dimension of (LST) databases for each group to compare the complexity of (LST) databases in different time intervals. This study revealed that the (LST) databases are mono-fractal, i.e. their spectrum displays a narrow width of scaling exponent, or in another word, a single global exponent is enough to characterize the fractal behavior and scaling properties of time series data. This important and novel finding helps to better understanding the structure of (LST) databases and consequently, better controlling the climate change caused by human activities. Although, because of limited number of databases, the vibration analysis using (PSD), fractal dimension and multi-fractal analysis could not classify (LST) databases for different seasons, however, they may be used as a comprehensive framework to further analysis, characterize and compare the fractal behavior of (LST) time series and extracting their complex structure. A future direction for this study would be utilizing larger number of databases from different climate, time scales and location to support our conclusion. plots, residuals do not follow a straight line very well and they deviate. As a result, the Land Surface Temperature (LST) data came from a Normal theoretical distribution or (LST) data sets are not Normally distributed (we can see bi-normality in figure (17)). Next, we plot Spread-Location plots of the Land Surface Temperature (LST) data (4) to check if we have equal variance or homoscedasticity in our database, i.e. if residuals are spreading equally or not. From figures ( 13)-( 16), the residuals are almost appeared to be randomly spread (do not spread widely). Finally, we plot the Residuals vs Leverage graphs of the Land Surface Temperature (LST) data (4) in figures ( 13)-( 16) to find influential subjects to determine the time series regression model. As we can easily see from the Residuals vs Leverage graphs, there is no obvious Cooks distance lines since all subjects are well located inside of the Cook's distance lines. Thus, we do not have any influential observation in (LST) databases. The diagnostic plots showed that the residuals of the linear model ( 1) are not normally distributed, therefore, the time series linear regression model may be unreliable to represent the Land Surface Temperature (LST) data (4). It is true that there exist many other different techniques such as R 2 , slope coefficients, and p-values to test the regression results, but, they may not be able to check all aspects of a perfect model.
However, the residuals can be considered as indicator to characterize if a model works perfectly or not. Moreover, residuals can well represent the unrevealed patterns in the data and need to be considered as a good platform to characterize all possible dynamics in different databases. In this study, we applied the time series linear regression to find the linear relationship between the (LST) variable Lt and regressor matrix Ht in model (1) using the Land Surface Temperature (LST) data in figure (4). In figure (12), we can see the time series linear regression analysis of the seasonal (LST) of Riley County during February, April, July and October (2021) obtained via Google Earth Engine in figure (4). Here, we analyze the time series linear regression in figure (12) to test if the linear model works well for the Land Surface Temperature (LST) data or not. At first, we check that if residuals have non-linear patterns or not. We have plotted the Residuals vs Fitted graph for February, April, July and October (2021) (4) in figures ( 13)-( 16) separately. According to the Residuals vs Fitted results, we have almost equally spreading residuals around the horizontal line (no distinct patterns), that means there is no non-linear relationships in any of time series regression models. Secondly, we are interested to check if residuals are normally distributed. Therefore, we plot Normal quantile-quantile (Q-Q) graphs of Land Surface Temperature (LST) data (4) in figures ( 13)-( 16) separately. In fact, Q-Q plots are scatter plots which are produced when we plot a set quantile against another one. As we can from the Q-Q In this study, we applied the time series linear regression to find the linear relationship between the (LST) variable Lt and regressor matrix Ht in model (1) using the Land Surface Temperature (LST) data in figure (4). In figure (12), we can see the time series linear regression analysis of the seasonal (LST) of Riley County during February, April, July and October (2021) obtained via Google Earth Engine in figure (4). Here, we analyze the time series linear regression in figure (12) to test if the linear model works well for the Land Surface Temperature (LST) data or not. At first, we check that if residuals have non-linear patterns or not. We have plotted the Residuals vs Fitted graph for February, April, July and October (2021) (4) in figures ( 13)-( 16) separately. According to the Residuals vs Fitted results, we have almost equally spreading residuals around the horizontal line (no distinct patterns), that means there is no non-linear relationships in any of time series regression models. Secondly, we are interested to check if residuals are normally distributed. Therefore, we plot Normal quantile-quantile (Q-Q) graphs of Land Surface Temperature (LST) data (4) in figures ( 13)-( 16) separately. In fact, Q-Q plots are scatter plots which are produced when we plot a set quantile against another one. As we can from the Q-Q