Spatial Distribution of Multi-Fractal Scaling Behaviours of Atmospheric XCO2 Concentration Time Series during 2010–2018 over China

Exploring the spatial distribution of the multi-fractal scaling behaviours in atmospheric CO2 concentration time series is useful for understanding the dynamic mechanisms of carbon emission and absorption. In this work, we utilise a well-established multi-fractal detrended fluctuation analysis to examine the multi-fractal scaling behaviour of a column-averaged dry-air mole fraction of carbon dioxide (XCO2) concentration time series over China, and portray the spatial distribution of the multi-fractal scaling behaviour. As XCO2 data values from the Greenhouse Gases Observing Satellite (GOSAT) are insufficient, a spatio-temporal thin plate spline interpolation method is applied. The results show that XCO2 concentration records over almost all of China exhibit a multi-fractal nature. Two types of multi-fractal sources are detected. One is long-range correlations, and the other is both long-range correlations and a broad probability density function; these are mainly distributed in southern and northern China, respectively. The atmospheric temperature and carbon emission/absorption are two possible external factors influencing the multi-fractality of the atmospheric XCO2 concentration. Highlight: (1) An XCO2 concentration interpolation is conducted using a spatio-temporal thin plate spline method. (2) The spatial distribution of the multi-fractality of XCO2 concentration over China is shown. (3) Multi-fractal sources and two external factors affecting multi-fractality are analysed.


Introduction
In recent decades, fossil fuel combustion and industrial processes arising from the intensification of anthropogenic activities have caused significant increases in the global atmospheric carbon dioxide (CO 2 ) concentration [1]. The excessively high atmospheric CO 2 concentration has led to a series of ecological and environmental problems, such as global warming, rising sea levels, droughts and floods, degradation of the soil ecosystem, and reductions in grain production [2][3][4]. Under the background of climatic change and global governance challenges, CO 2 , as the most important and longest-lived greenhouse gas in the atmosphere, has attracted increasing attention from researchers and policymakers. Analyses of the spatial distributions of the multi-fractal scaling behaviours, including small or large fluctuations, multi-fractal strength and multi-fractal sources, of the atmospheric CO 2 concentration records are helpful for understanding the dynamic mechanisms of CO 2 concentration variations owing to carbon emission/absorption and climatic change.
In the past several years, an ever-growing number of studies have focused on the multi-fractal scaling behaviours of meteorological and air quality parameters in the field of atmospheric sciences using the well-established multi-fractal detrended fluctuation Entropy 2022, 24, 817 2 of 14 analysis (MFDFA) [5,6]. The MFDFA has been applied into time series concerning atmospheric temperature [7][8][9][10][11], precipitation [12,13], relative humidity [14], particulate matter concentrations [15], and many other measurements. A few investigations have focused on spatial distributions of the multi-fractal scaling behaviours for these meteorological elements [10,16,17]. These investigations, especially those concerning the spatial heterogeneity of meteorological elements, are more helpful for understanding the corresponding dynamic mechanisms.
However, to present, in the context of the analysis of the complexity and non-linearity of signals, the multi-fractal scaling behaviours of CO 2 mainly concern two aspects: the carbon trade, and the atmospheric CO 2 concentration record (with various time scales including hour, day, and month). In the literature, there are numerous investigations focusing on the scaling behaviours of the carbon trade [18][19][20][21][22][23][24]; however, relatively fewer investigations have focused on the atmospheric CO 2 concentration [25][26][27][28], particularly from the perspective of spatial heterogeneity. Evidently, our knowledge of the scaling behaviours of atmospheric CO 2 concentrations is limited. For instance, some researchers believed that atmospheric CO 2 concentration records were characteristic of a long-term memory and obeyed a power law [25]; in contrast, others argued that thermospheric CO 2 power did not exhibit power-law behaviours [28]. Their standpoints contradict each other.
As the direct sampling of greenhouse gas concentrations, especially in the upper atmosphere, requires considerable effort and massive costs [29], current ground-based measurement stations are sparse and uneven [29]. Spatio-temporal CO 2 concentration records are not plentiful. Our knowledge of the scaling behaviours of CO 2 concentration is further limited by the insufficiency of spatio-temporal CO 2 concentration records.
Nonetheless, Japan launched the world's first greenhouse gas observation satellite, the Greenhouse Gases Observing Satellite (GOSAT) in 2009 [30,31]. This satellite has offered a potential way to complement existing measurements of CO 2 concentrations through space-based remote sensing. Since CO 2 concentration values from GOSAT have been released, they have been practically applied in many scientific investigation works [32][33][34][35][36][37]. Although space-based observations can obtain CO 2 concentrations for global coverage and have a high measurement density, these satellite observations have gaps and are irregularly positioned owing to certain limitations including clouds, the observational mode of the satellite [30], and inversion algorithms. It is necessary to fill these gaps, and geospatial statistics might be an effective tool for this task [38]. Ref. [39] incorporated temporal variability to improve a Kriging geo-statistical analysis [40] of satellite-observed CO 2 concentration in China, and presented a reliable interpolation of the XCO 2 concentration. However, in the spatio-temporal Kriging interpolation [40], it was very difficult to select a reliable and reasonable semi-variance model and its parameter [40]. Inappropriate selections of models and parameters can lead to significantly improper results and low efficiencies in interpolation.
In this study, to overcome the lack of knowledge regarding the multi-fractal scaling behaviour, we investigate the spatial distribution of the multi-fractal scaling behaviour of an atmospheric XCO 2 concentration time series over China from 2010−2018, and get an insight into the dynamic mechanisms of CO 2 concentration changes. China covers a large geographical area, with wide climatic variety, complicate topography and plentiful vegetation coverage types. Within this context, it would be interesting and ideal to examine the spatial distribution of the multi-fractal properties over the whole of China. We propose an improved spatio-temporal interpolation approach to implement spatio-temporal interpolation of atmospheric XCO 2 concentration, and use MFDFA to estimate their multi-fractal scaling exponents. The results of this study are respected to contribute to the knowledge regarding the multi-fractal scaling behaviours, their spatial distribution, multi-fractal sources, and possible external influencing forces.

Study Area
The study region covers the Chinese mainland with latitude from 21 • to 54 • N and longitude from 70 • to 138 • E, divided into a total of 4189 grid points with a spatial resolution of 0.5 • × 0.5 • . Seas and islands in China are excluded, as the greenhouse gas tendency largely varies between land and sea. Taiwan, Hong Kong, and Macao are also excluded, owing to an absence of data sources. Figure 1 shows spatial pattern of major vegetation types (http://data.tpdc.ac.cn/zhhans/data/eac4f2cf-d527-4140-a35d-79992957f043/ (accessed on 10 June 2022)). The major land-cover types in China are forest, grassland, paddy fields and rain-fed croplands, which are mainly located in northeast China, northern China, and southern China, respectively. More specifically, a deciduous broadleaf forest is located in northeast China, and an evergreen broadleaf forest is distributed in southern China. Most of the land-cover types in the earth's terrestrial ecosystems are found in China [41].
cluded, owing to an absence of data sources. Figure  1 shows spatial pattern of major vegetation types (http://data.tpdc.ac.cn/zh-hans/data/eac4f2cf-d527-4140-a35d-79992957f043/ (accessed on 10 June 2022)). The major land-cover types in China are forest, grassland, paddy fields and rain-fed croplands, which are mainly located in northeast China, northern China, and southern China, respectively. More specifically, a deciduous broadleaf forest is located in northeast China, and an evergreen broadleaf forest is distributed in southern China. Most of the land-cover types in the earth's terrestrial ecosystems are found in China [41].
There are four typical grid points with different climate, land cover and vegetation types, marked by red coloured circles (A, B, C and D). The grid point A (109°15′ E, 42°15′ N) is located in Inner Mongolia, with low air temperature and little precipitation, where land is mainly covered by grassland or barren land like desert. The grid point B (116°45′ E, 23°15′ N) is placed in Guandong province, with high air temperature and plentiful precipitation, where major land cover is forestry. The grid point C (120°45′ E, 32°45′ N) is situated in Jiangsu province, with median air temperature and moderate precipitation. The double crop system is the main characteristic of this region. The grid point D (100°54′ E, 36°17′ N, 3816-m altitude) is WLG (ground observational station), with low air temperature and high altitude. WLG is the highest atmospheric background reference monitoring station in the world. As WLG station is far away from the built-up area, the CO2 concentration measured by an instrument is highly accurate.  is WLG (ground observational station), with low air temperature and high altitude. WLG is the highest atmospheric background reference monitoring station in the world. As WLG station is far away from the built-up area, the CO 2 concentration measured by an instrument is highly accurate. The Greenhouse Gases Observing Satellite was launched in 2009 by Japan, and circles around the earth in a sun-synchronous orbit at a 666-km altitude with the descending node at approximately 12:48 local time and a 3-day recurrence [30,31]. GOSAT was carefully and specifically designed to record two concentrations of greenhouse gases from space: the column-averaged dry-air mole fraction of carbon dioxide (XCO 2 ) and methane (XCH 4 ), respectively [30]. The GOSAT Level 2 data are the CO 2 dry air mixing ratios (XCO 2 ), defined by the ratio of the total number of CO 2 molecules to that of dry molecules, not only near the earth surface, but also in the entire vertical column to the top of the atmosphere [30].
We acquired and rearranged the monthly XCO 2 data from GOSAT from January 2010 to December 2018. A total of 10367 observations were available. Theoretically, a total of 452,412 (4189 grid points × 9 years × 12 months per year) observations should be available for all of the grid-points in the study region. Clearly, various observations in the spatial and temporal dimensions were missing. For instance, all of the CO 2 concentration values in January 2015 were absent. In an additional case, the CO 2 concentration values in May 2010 were available, but the ones in July 2010 were missing at the same location. Therefore, a spatio-temporal thin plate spline interpolation approach (STTPS) needed to be developed.

In-Situ CO 2 Concentration Values from Ground-Based Station
To examine the interpolation accuracy of XCO 2 concentration, the in-situ CO 2 concentration values were collected from the Waliguan meteorological ground-based observation station. The atmospheric environment and components of the WLG station were very stable, owing to its distance from built-up areas. Therefore, the precision and reliability of the in-situ observed values of the CO 2 concentration at WLG station are high. Notably, the CO 2 concentration values for April 2011, November 2011, December 2011, and August 2014 were absent from the WLG during 2010−2018.

Spatio-Temporal Thin Plate Spline Interpolation
As the number of observations of spatio-temporal XCO 2 concentration is finite, it is difficult to use these XCO 2 concentration observations to realize the valuable analysis. The spatio-temporal interpolation, thus, is a simple and economical approach for obtaining sufficient observations of the spatio-temporal XCO 2 concentration.
Let us suppose that in the study region, there is a set of CO 2 data values in space-based grid points. To obtain a continuous surface consisting of uniformly distributed grid points with a spatial resolution of 0.5 • × 0.5 • , the interpolation is conducted based on a radial basis function (RBF) [42,43]. The RBF interpolation (a distance-based function) is written as Equation (1), as follows: In the above, d j is a Euclidean distance between an interpolated point and the jth known point, and respective x and y are the latitude and longitude coordinates of the interpolated point. Among the distance-based functions, the thin plate spline function is the most fundamental type of interpolation one. It is considered as a bending thin plate, and can pass through all known points with minimal bending energy. Hence, in Equation (1), there are a total of n + 3 interpolation coefficients, namely, A j (j = 1, 2, 3 . . . , n), b 1, b 2 , and b 3 that need to be estimated. This is called as the spatial thin plate spline interpolation method (abbreviated as STPS). To overcome singularity when estimating the interpolation coefficients, three additional constraint conditions derived from the spline function properties should be incorporated as follows: When Equations (1) and (2) are combined, all the interpolation coefficients can be estimated easily. Next, temporal variability is incorporated to improve STPS interpolation of the satellite-observed CO 2 measurements. Equation (1) is thus rewritten into Equation (3) as follows: Here, ϕ(x, y, t) is an interpolated CO 2 concentration, t is the time, and d 2 j log d j is a RBF with the minimum spatio-temporal surface curvature. d j is a spatio-temporal distance between the interpolated point and the jth known point. It is denoted by d 2 and τ is a spatio-temporal ratio. There are the interpolation coefficients of n + 4 (A j (j = 1, 2, 3 . . . , n), b 1 , b 2 , b 3 , and b 4 ). Four additional constraint conditions analogous to those in Equation (2) need to be incorporated before estimating the interpolation coefficients. This is called as the spatio-temporal thin plate spline interpolation method.
The optimal τ, denoted as τ*, is determined by experiments as follows. First, 90% of the known data points are randomly chosen from the XCO 2 concentration dataset as the interpolation control points for estimating the interpolation coefficients, and the remaining 10% data points are used to examine interpolation accuracy. Second, a different value of τ is set to get a different interpolation variance. A plot of the interpolation variance against τ is drawn. When the interpolation variance is the smallest, the corresponding τ is then set as τ*. To obtain τ* as fast as possible, the dichotomy is generally adopted.

Multi-Fractal Detrended Fluctuation Analysis
The generalised Hurst exponent (scaling exponent) is important for revealing the long-term correlation in time series [44,45]. A detrended fluctuation analysis is a powerful and effective tool for estimating the generalised Hurst exponent of non-stationary time series. Here, we briefly recall the detrended fluctuation analysis (DFA) [44] and MFDFA [5].
Assume that there is one time series {x(k)}, k = 1, 2, . . . , N, where N is the length of this series. First, one must calculate the profile, as follows: In the above, x denotes the averaging over the whole time series {x(k)}. Next, one needs to divide the profile Z(i) into N s = int (N/s) non-overlapping segments of equal length s. As the series' length N is not always a multiple of the given time scale s, a short part may remain at the end of each profile. To avoid discarding this part of the series, the same procedure is repeated starting from the opposite end of each profile, and 2N s segments are thus obtained together. The detrended covariance is determined as follows: Here, z λ (j) is the fitting polynomial in box λ. Next, the fluctuation function is defined in Equation (6) [44]. Equation (6) can be rewritten into Equation (7), as follows: The varying moment q generally takes any real value. Similarly to the above, the linear relationship on a log-log plot exhibits the existence of a power-law F q (s) ∝ s h(q) . If the time series is mono-fractal, h(q) is independent of q; in contrast, if the time series is multifractal, h(q) depends on q. For q > 0, h(q) describes the behaviour of the boxes with large fluctuations, whereas for q < 0, h(q) depicts the behaviour of the boxes with small fluctuations. Via a Legendre transform, we can obtain singularity strength and fractal dimension in Equations (8) and (9), as follows: In the above, τ(q) = qh (q) −1, the quantity α is singularity strength (or Hölder exponent), whereas f (α) is a singularity spectrum for expressing the dimensionality of the subset of the time series [5]. A large ∆α indicates that the multi-fractality is strong, whereas a small ∆α suggests that the multi-fractality is weak. There are two types of multi-fractal sources: (i) a broad probability density function, and (ii) long-range correlations [5]. To distinguish multi-fractality sources, we might use a shuffle series and a surrogate series (or phase-randomised series). If the shuffle series has any multi-fractality, it would be owing to the broad probability density function; whereas if the surrogate series exhibits multifractality, it would be owing to the long-range correlations. In nature, the multi-fractality is sometimes driven by one of the sources or by both simultaneously.

Accuracy Evaluation of the Interpolated Monthly XCO 2 Concentration
The STTPS is applied to interpolate the monthly XCO 2 concentrations from January 2010 to December 2018 over China. To examine the interpolation accuracy of XCO 2 concentration, the average over the square of difference between the interpolated value and observed value is denoted by σ 2 . The smaller the value of σ 2 , the higher the interpolation accuracy is.
The two curved lines of the in-situ observed values and interpolated values at WLG (the ground-based station) during the identical period are shown in Figure 2a, where the black solid curved line is the observed in-situ value, and the dotted line is the interpolated one. Their σ 2 is 7.88. As a whole, the trends of the two curved lines are in basic agreement, reflecting the identical seasonal fluctuations and inter-annual growth of the XCO 2 concentration.
To confirm that STTPS has the advantage over STPS interpolation, the STPS interpolation method is also applied for the same data sources. Figure 2b exhibits the comparison between the observed values and interpolated values at WLG. Their σ 2 is 9.95. Apparently, the value of σ 2 when using STTPS is lower than that when using STPS approach. Furthermore, the peaks in the interpolated values when using STPS approach and those To confirm that STTPS has the advantage over STPS interpolation, the STPS interpolation method is also applied for the same data sources. Figure 2b exhibits the comparison between the observed values and interpolated values at WLG. Their σ 2 is 9.95. Apparently, the value of σ 2 when using STTPS is lower than that when using STPS approach. Furthermore, the peaks in the interpolated values when using STPS approach and those in the observed values do not appear simultaneously in Figure 2b, whereas the peaks in the interpolated values using STTPS and those in the observed values occur almost simultaneously in Figure 2a. Therefore, the interpolated values using STTPS have rather high interpolation accuracy. To further evaluate the interpolation accuracy, a ten-fold cross validation method was employed. In particular, 90% of the satellite-based CO2 data points were randomly selected as interpolation control points, and the remaining 10% were used as interpolation validation points. After the optimal τ * was determined, STTPS was conducted. The identical interpolation process was repeated 100 times, and the average values of the evaluation indicators, such as the mean absolute error (MAE), mean square error (MSE), and root mean square interpolation error (RMSE), were calculated. If the MAE, MSE and RMSE are close to zero, it means that the interpolated values are approximate to the observed values, i.e., the interpolation accuracy is high. Table 1 shows the respective averages of MAE, MSE, and RMSE using STTPS. Our results were compared with those in Ref. [39]. It is intuitively inferred that interpolation values are consistent from the perspective of a thematic map, and that the values σ 2 are roughly equivalent to each other. The interpolation accuracy is good.  To further evaluate the interpolation accuracy, a ten-fold cross validation method was employed. In particular, 90% of the satellite-based CO 2 data points were randomly selected as interpolation control points, and the remaining 10% were used as interpolation validation points. After the optimal τ* was determined, STTPS was conducted. The identical interpolation process was repeated 100 times, and the average values of the evaluation indicators, such as the mean absolute error (MAE), mean square error (MSE), and root mean square interpolation error (RMSE), were calculated. If the MAE, MSE and RMSE are close to zero, it means that the interpolated values are approximate to the observed values, i.e., the interpolation accuracy is high. Table 1 shows the respective averages of MAE, MSE, and RMSE using STTPS. Our results were compared with those in Ref. [39]. It is intuitively inferred that interpolation values are consistent from the perspective of a thematic map, and that the values σ 2 are roughly equivalent to each other. The interpolation accuracy is good. There are a total of 4189 grid points in the study region (See Figure 1). We analysed the multi-fractal behaviour of the month-average XCO 2 concentration time series for every grid point. For convenience of discussion, we first analyzed multi-fractal behaviours of four typical grid points.
The temporal evolution and profile of the original month-average of XCO 2 concentration time series at the WLG ground observational station from 2010 to 2018 are shown in Figure 2a. Apparently, the trends and approximate year-periodical fluctuations are easily observed. To reduce the multi-fractality uncertainties caused by the year-periodical fluctuations, the global year-periodical fluctuation should be removed by using a fitted model, as follows: (10) In the above, x(t) is a month-average XCO 2 concentration, β is a slope, C is an intercept, res(t) is a remainder component, and period 1 -period 3 are the approximate annual period, half-annual period, and seasonal period, respectively. Ref. [7] recommended using a seasonal-trend decomposition procedure based on Loess (STL) to decompose a time series [46]. We compared the periodicity strength of the remainder component in Equation (10) with that from STL method using the wavelet analytical technology. As a result, the periodical strength of the remainder component in Equation (10) is too small to be seen. The fitted model in Equation (10) is more appropriate than the STL method. DFA3 is recommended to estimate generalised Hurst exponent h(q) in MFDFA. Figure 3 shows the generalised Hurst exponents of the original series, shuffle series, and surrogate series against q for four typical grid points A, B, C, and D, respectively. Herein, the panels (A1−D1) express respective h(q) of four original series with q, panels (A2−D2) are respective h suf (q) of four shuffle series against q, and panels (A3−D3) are respective h sur (q) of four surrogate series versus q. The q ranges from −0.5 to 2.5, as the original series' length is somewhat short. As seen, the exponents h(q) of the original series generally decrease with the increasing q. The respective value ∆h = max(h(q)) − min(h(q)) in the panel (A1 and D1) is relatively small, indicating that their multi-fractality is weak in the grid point A and D; whereas the respective value ∆h in panel (B1 and C1) is comparatively large, revealing that their multi-fractality is strong in the grid point B and C.     panels (b and c) show a strong multi-fractality, whereas the panels (a and d) display a weak multi-fractality. This finding is in agreement with Figure 3. In addition, the two singular spectra in panels (b and c) are asymmetrical, indicating that their multi-fractality dominantly stem from a large fluctuation of the original series.
To distinguish between the multi-fractal sources for the above four typical grid points, we first shuffled the original series and then estimated h suf (q) with q from −0.5 to 2.5. Next, we repeated the above procedure 1000 times, and ultimately averaged h suf (q) and calculated variance of h suf (q)(denoted as σ h (q)). Consequently, the maximum variance σ h (q) for the grid points (A and D) is 0.15. The panels (A2-D2) in Figure 3 show the averaged h suf (q) versus q. ∆h suf = max(h suf (q)) − min(h suf (q)) is defined. We judged that if ∆h suf < 0.15, the shuffle series is mono-fractal, indicating Entropy 2022, 24, 817 9 of 14 that its multi-fractality is owing to long-range correlations; otherwise the multi-fractal source is a broad probability density function. This criterion is a bit crude, as the original series length is finite. Figure 3. Generalised Hurst exponents of the original series, shuffle series, and surrogate se against q for four typical grid points A, B, C, and D, respectively. Panels (A1A3) for grid poin panels (B1B3) for grid point B, and so on. To distinguish between the multi-fractal sources for the above four typical g points, we first shuffled the original series and then estimated h suf (q) with q from −0.5 2.5. Next, we repeated the above procedure 1000 times, and ultimately averaged h su and calculated variance of h suf (q)(denoted as σh(q)). Consequently, the maximum varia σh(q) for the grid points (A and D) is 0.15. The panels (A2-D2) in Figure 3 show the eraged h suf (q) versus q. Δh suf = max(h suf (q)) − min(h suf (q)) is defined. We judged that if Δ < 0.15, the shuffle series is mono-fractal, indicating that its multi-fractality is owing long-range correlations; otherwise the multi-fractal source is a broad probability dens function. This criterion is a bit crude, as the original series length is finite. Subsequently, similarly to the shuffle procedure, we repeatedly generated the surrogate series 1000 times by using phase-randomization, then estimated and ultimately averaged h sur (q). The panels (A3-D3) in Figure 3 show the averaged h sur (q) versus q. We judged that if ∆h sur < 0.15, the surrogate series is mono-fractal, and the multi-fractal source is a broad probability density function; otherwise, the multi-fractal source is long-range correlations.
As shown in the panels B2, B3, C2 and C3 in Figure 3, we discerned that the multi-fractal source originates from long-range correlations, and not from a broad probability density function. We determine whether the remainder components of the two original series in grid points B and C are characteristic of a long-tail distribution using the Statistical Package for Social Sciences. As a result, the two remainder components do not pass the statistical test for the long-tail distribution under a confidence interval of 95%.

Spatial Distribution of Multi-Fractal Scaling Behaviour
The spatial distribution of the scaling exponents H = h(q = 2) for all of the original series is shown in Figure 5a. The values of h(q = 2) are greater than 0.6 in most regions in China, indicating that XCO 2 concentration records in most regions show a persistent auto-correlation. This finding suggests that the XCO 2 -concentration fluctuation continues to increase from 2010 to 2018. However, in a few regions like the central part of Shanxi province, the values of h(q = 2) range from 0.44 to 0.55, indicating that the XCO 2 -concentration series are uncorrelated. Figure 5b shows the spatial distribution of multi-fractal strength ∆α of the original series. When ∆α < 0.3, the original series are mono-fractal otherwise are multi-fractal [10]. As a whole, the multifractal behaviour in western China is stronger than that in eastern China. Within eastern China, the multi-fractal behaviour in the southern part is stronger than that in northern part, whereas in the upstream areas between the Yangtze River and Yellow River and along the Yangtze River Basin, the multi-fractal behaviour is weak. We notice that there are few grid points with ∆α ≈ 0, i.e., less than 1% of the total grid points.

Spatial Distribution of Multi-Fractal Scaling Behaviour
The spatial distribution of the scaling exponents H = h(q = 2) for all of the original series is shown in Figure 5a. The values of h(q = 2) are greater than 0.6 in most regions in China, indicating that XCO2 concentration records in most regions show a persistent auto-correlation. This finding suggests that the XCO2-concentration fluctuation continues to increase from 2010 to 2018. However, in a few regions like the central part of Shanxi province, the values of h(q = 2) range from 0.44 to 0.55, indicating that the XCO2-concentration series are uncorrelated.    The spatial distribution of the multi-fractal sources is shown in Figure 6. There are mainly two types of multi-fractal sources: one is long-range correlations, and the other is both long correlations and a broad probability density function. The former is mainly distributed in southern, middle and northeast China, and accounts for 51% of the grid points. The approximately triangle-shaped pattern is shown in southern and middle China (light blue region in Figure 6). The latter is located in northern and middle China with an inverted-triangle-shaped pattern (yellow region part in Figure 6), and accounts for the remaining 44% of all of the grid points.
In summary, the scaling behaviours are a positive and persistent auto-correlation in most regions. The multi-fractal strength of the original XCO 2 concentration series is different, as it is strong in western China and weak in eastern China. There are two multi-fractal sources: one is long-range correlations and the other is both long-range correlations and a broad probability density function [5].
both long correlations and a broad probability density function. The former is mainly distributed in southern, middle and northeast China, and accounts for 51% of the grid points. The approximately triangle-shaped pattern is shown in southern and middle China (light blue region in Figure 6). The latter is located in northern and middle China with an inverted-triangle-shaped pattern (yellow region part in Figure 6), and accounts for the remaining 44% of all of the grid points. Figure 6. Spatial distribution of the multi-fractal sources. P-P is a broad probability density function, L-P is both long-range correlations and a broad probability density function, and L-L is long-range correlations.
In summary, the scaling behaviours are a positive and persistent auto-correlation in most regions. The multi-fractal strength of the original XCO2 concentration series is different, as it is strong in western China and weak in eastern China. There are two mul- Figure 6. Spatial distribution of the multi-fractal sources. P-P is a broad probability density function, L-P is both long-range correlations and a broad probability density function, and L-L is long-range correlations.

Discussions
The multi-fractal strength in China is different, i.e., strong in western China and weak in eastern China. The CO 2 concentration in the atmosphere is usually controlled by two physical mechanisms: climatic changes like atmospheric temperature, and carbon emission/absorption. Owing to two physical mechanisms, the dynamic mechanisms of the CO 2 concentration seem more complicated than a single physical mechanism.
Refs. [9,10] reported that the strong multi-fractality of atmospheric temperature was mainly located in southern and northern China, while the weak multi-fractality was along the Yangtze River Basin. By comparison, the spatial distribution of the multi-fractality of XCO 2 concentration in Figure 5b somewhat resembles that of the atmospheric temperature [9,10]. This finding is easily understood. An increase in atmospheric temperature will cause CO 2 molecules to leave open systems out. The increasing atmospheric temperature can thereby lower the CO 2 concentration. We can, hence, preliminarily infer that the atmospheric temperature, as an external force, has an effect on the multi-fractality of the XCO 2 concentration, i.e., has an influence on small fluctuations and large fluctuations of XCO 2 concentration. Figure 5c shows that the small fluctuations of XCO 2 concentration in southern China are weaker than in northern China. We speculated that vegetation absorption is one of external forces causing the small fluctuations of XCO 2 concentration. The absorption magnitude and absorption time length have strong influences on small fluctuations. We noticed that vegetation coverage in eastern China is high, whereas it is very low in western China owing to the desert. Although most south areas of the Yangtze River and the northeast region (eastern region in Inner Monglia, Heilongjiang and Jilin) in China are covered by forestry, growing season of forestry is different. The growing season in south areas of the Yangtze River is usually from April to October even November in some places, but in northeast region, is from May to September. Clearly, the length of absorption time of CO 2 in south areas of the Yangtze River is longer than northeast region. The number of small fluctuations in the original series for XCO 2 concentration in south areas of the Yangtze River thereby is more than northeast region, indicating that the values α max in south areas of the Yangtze River are greater than northeast region. Meanwhile, in the upstream areas between the Yangtze River and Yellow River and along the Yangtze River Basin, double cropping system has influences on fluctuations of CO 2 concentration to some extent. The wheat harvest is frequently completed at the end of May every year in this region, and then corn or soybeans are planted immediately in early June. Two absorption peaks of CO 2 concentration often appearing in May and July, respectively, are clearly discerned. Hence, compared with south areas of the Yangtze River, the number of small fluctuations is less, and the values α max is small.
The spatial differentiation of the vegetation cover has a strong effect on CO 2 concentration. Besides that, human-induced disturbances including industry production, vehicle emissions, and nature-induced variations like local forest fires, biomass burning, and cultivated land destruction, will affect the variation in CO 2 concentration owing to emissions. Thereby, we speculate that it is possible that two external forces affect the multi-fractality of the CO 2 concentration: one is climatic change like atmospheric temperature, and the other is carbon emission/absorption. However, the specific physical mechanisms as to how these two external forces affect the multi-fractal sources of the CO 2 concentration merit further investigation.

Conclusions
In this work, we investigated the spatial distributions of multi-fractal scaling behaviours of an interpolated XCO 2 concentration time series from 2010-2018 at a total of 4189 grid points by applying MFDFA method. The main results can be summarised as follows.

(I)
We improved a spatio-temporal thin plate spline interpolation approach, and conducted interpolation of the monthly XCO2 concentrations over China from 2010−2018 based on GOSAT observations of XCO2. The interpolation accuracy of spatio-temporal thin plate spline interpolation approach was higher than a spatial thin plate spline interpolation one. The interpolated XCO2 concentration is highly accurate and is useful in analyzing multi-fractal scaling behaviours. (II) We found that the scaling behaviours of XCO2 concentration show a positive and persistent autocorrelation in most regions. The scaling behaviours of CO2 did not always obey power laws. (III) The multi-fractal strength of XCO2 concentration is different, i.e., strong in western China and weak in eastern China. There are two types of multi-fractal sources: one is long-range correlations, and the other is both long-range correlations and a broad probability density function. Two types are mainly distributed in southern and middle China with a triangleshaped pattern and in northern China with an inverted-triangle-shaped pattern, respectively. Two external forces are likely to have influences on multi-fractality: the climatic change like atmospheric temperature and the carbon emission/absorption.
There is a limitation in this study, that is, the original XCO 2 concentration time series length is relatively short. Because the spatial and temporal resolutions of the XCO 2 concentration are low, the data sources directly available in the investigation are insufficient. If we use fewer observational data sources to conduct the spatio-temporal interpolation, the uncertainty in the interpolated data would greatly increase. Thus, mining suitable data sources of long CO 2 concentration series with high quality (direct observational values instead of interpolation ones) is another important task in the future.