Modeling Soil CO₂ Efflux in a Subtropical Forest by Combining Fused Remote Sensing Images with Linear Mixed Effect Models

Abstract

Monitoring tropical and subtropical forest soil CO₂ emission efflux (FSCO₂) is crucial for understanding the global carbon cycle and terrestrial ecosystem respiration. In this study, we addressed the challenge of low spatiotemporal resolution in FSCO₂ monitoring by combining data fusion and model methods to improve the accuracy of quantitative inversion. We used time series Landsat 8 LST and MODIS LST fusion images and a linear mixed effect model to estimate FSCO₂ at watershed scale. Our results show that modeling without random factors, and the use of Fusion LST as the fixed predictor, resulted in 47% (marginal R² = 0.47) of FSCO₂ variability in the Monthly random effect model, while it only accounted for 19% of FSCO₂ variability in the Daily random effect model and 7% in the Seasonally random effect model. However, the inclusion of random effects in the model’s parameterization improved the performance of both models. The Monthly random effect model that performed optimally had an explanation rate of 55.3% (conditional R² = 0.55 and t value > 1.9) for FSCO₂ variability and yielded the smallest deviation from observed FSCO₂. Our study highlights the importance of incorporating random effects and using Fusion LST as a fixed predictor to improve the accuracy of FSCO₂ monitoring in tropical and subtropical forests.

1. Introduction

Forest soil CO₂ emission (FSCO₂) serves as a crucial conduit in the global carbon cycle and exerts a significant impact on the carbon budget of terrestrial ecosystems. In recent decades, the continuous increase in greenhouse gases in the atmosphere has changed forest growth and productivity, resulting in the acceleration of carbon dioxide emissions from soil to the atmosphere [1]. Therefore, timely and accurate observation of FSCO₂ at large spatial scales has significant scientific value for deeply revealing the mechanism of the global carbon cycle and for precisely predicting future climate change.

There are many FSCO₂ monitoring methods, among which the Box method is widely used [2]. Although the automatically closed gas chamber can effectively obtain the change of FSCO₂ in time, most gas chambers still need to fix the spatial difference of FSCO₂ monitoring at the regional or landscape scale. In addition, due to the lack of long-time-series field observation of FSCO₂, earth observation images cannot directly monitor the temporal and spatial variation of FSCO₂, which limits the direct prediction of FSCO₂ through remote sensing products [3]. Therefore, previous studies estimated FSCO₂ indirectly through biotic and abiotic factors, among which soil temperature is the key environmental factor to control FSCO₂ [4]. These studies mainly investigated interconnections between soil temperature and carbon release from the land surface from different aspects [5,6,7].

Remote-sensing-based land surface temperature (LST) measurement has great potential for obtaining the spatial distribution characteristics of FSCO₂ at the pixel level [1,8,9,10]. The MODIS LST product is one of the successful examples. It is currently widely used in monitoring large-scale FSCO₂ around the world [11]. The split window algorithm, used for inverting this product, can effectively compensate for atmospheric attenuation. However, the limitation in spatial resolution of MODIS LST restricts its applicability at small, regional scales. Landsat 8 addresses the shortcomings of MODIS LST products and provides an unprecedented opportunity to improve the accuracy of LST measurement via remote sensing inversion. Several case studies have used the split window algorithm to obtain LST from Landsat 8 [12,13,14,15]. However, a recent report shows that the thermal band 11 of Landsat 8 is subject to relatively high levels of stray light interference, which induces caution in using a split window algorithm [16]. Therefore, the data fusion algorithm has been used to estimate LST from single-channel TIRS images [17,18,19,20]. Although the quality of LST data with different temporal and spatial resolutions has been improved greatly, the mapping of FSCO₂ based on the improved MODIS and Landsat 8 LST remote sensing fusion data is still scarce.

The aim of this study is to address the following scientific problems in estimating the FSCO₂ at high spatiotemporal scale: first, how to retrieve a high-quality spatiotemporal data set of remote sensing-based FSCO₂ observation variables; second, how to construct an efficient estimation model using multisource satellite-based data combined with an FSCO₂ model. The solution to this problem is valuable for accurately estimating the global carbon budget, mitigating global warming, and accurately predicting future climate change.

2. Materials and Methods

2.1. Study Basin

The study area (23.67–23.96°N, 114.03–113.75°E) is a headwater catchment of the Liuxihe (LXH) River basin, named the LXH Reservoir watershed, and located in the northern part of Guangdong Province of China [21] (Figure 1). It is characterized by a typical subtropical monsoon climate with long-term mean annual temperature of 20.3 °C and precipitation of 2100 mm [22]. The rainfall is concentrated in the period from March to September (about 80%). The rest of the months are part of the dry period in this basin, mainly consisting of autumn and winter. The catchment covers 456.7 km², and the altitude ranges from about 150 m to more than 1140 m. As one of the major head water reservoir basins in southern China, it supplies drinking water to the city of Guangzhou and thus is an important ecological barrier in the north of Guangzhou [23]. Vegetation cover, mainly forests, accounts for 80% (365.4 km²) of the catchment area. The main vegetation types are broad-leaf and needle-leaf evergreen forests. We selected Chenhedong (CHD) Nature Reserve and GuoYuan (GY) in the LXH National Forest Park as sampling points and conducted 10 independent measurements of the FSCO₂ at each.

2.2. Data Processing

2.2.1. Earth Observation Data Sets and Preprocessing

The remote sensing data employed in this study consist of Landsat 8 images and Land Surface Temperature (LST) products derived from daily Moderate Resolution Imaging Spectroradiometer (MODIS) images. The MODIS LST product (MOD11A2) was obtained from the Oak Ridge National Laboratory Distributed Active Archive Center (DAAC) (http://MODIS.ornl.gov/cgi-bin/MODIS) while the Landsat 8 images were downloaded from the US Geological Survey (USGS) website for free (http://Landsat.usgs.gov). A total of 70 MODIS LST products and 7 Landsat images (

Table 1. Metadata for Landsat 8 images.

Date	Cloud Cover (%)	Azimuth/Zenith Angles of the Sun	ML(10)	AL(10)	K1	K2
11 September 2019	0	135.06/58.8	3.34 × 10−4	0.1	774.89	1321.08
27 September 2019	5.9	139.15/57.0	3.34 × 10−4	0.1	774.89	1321.08
29 October 2019	1.5	151.85/47.8	3.34 × 10−4	0.1	774.89	1321.08
14 November 2019	0.3	154.61/43.4	3.34 × 10−4	0.1	774.89	1321.08
30 November 2019	52.3	155.31/39.9	3.34 × 10−4	0.1	774.89	1321.08
17 January 2020	91.3	149.15/38.2	3.34 × 10−4	0.1	774.89	1321.08
18 February 2020	0	141.45/45.3	3.34 × 10−4	0.1	774.89	1321.08

M_L(10) is the multiplicative rescaling factor; A_L(10) is the additive rescaling factor; K₁ and K₂ are thermal constants; Cloud cover represents the cloud coverage of a scene of Landsat 8 image.

) with low cloud coverage in the study region (>1%) were used.

We collected MODIS LST products of day and night spanning from 28 September 2019 to 4 January 2020. Filtering and clipping of images were complemented before usage. Landsat 8 images without cloud cover in the study area were collected for the study period (row 122 and path 43). In order to counteract the influence of stray light on the thermal band 11 of Landsat 8, this study utilized thermal band 10 (with a spatial resolution of 100 m) for estimating LST.

2.2.2. Field-Observed Data and Preprocessing

The field measurements of FSCO₂ were conducted using a closed chamber in each study site (GY and CHD) during the study period (Figure 1). FSCO₂ at an interval of two times per month for each site was measured using an automatic cavity ring-down spectrophotometer [22]. There was no special setting for the instrument in this study, and the default company settings are kept for reference [24].

Furthermore, in situ measurements of soil temperature, soil moisture, air temperature, and precipitation were also conducted. Specifically, soil temperature (ST) and moisture (SM) were measured at soil depths of 5–10 cm in the vicinity of the chamber sampling site using a heat dissipation probe (TDP, Dynamax Inc., Houston, TX, USA) at three different locations [22]. The mean values of these variables were obtained by averaging repeated measurements at each site.

2.3. Methods

Figure 2 illustrates the specific technical methods used in this study. To improve the precision and accuracy of monitoring the carbon release flux from forest soils in the study area, high-quality spatiotemporal remote sensing fusion images were obtained, and the LXH headwater watershed in the humid region of South China was selected as the research area. High temporal and spatial remote sensing forest soil CO₂ emission flux parameters were retrieved, and the new fusion algorithm and the CO₂ emission flux estimation model were effectively coupled to construct the multi-state variable optimization model.

Finally, the field measurements of FSCO₂ and high-quality spatiotemporal remote sensing parameter sets of carbon flux were used for accuracy assessments of the forest soil CO₂ emission in different models, which improved the accuracy and efficiency of regional CO₂ emission flux estimation.

2.3.1. Landsat 8 LST Calculation

In order to analyze the data effectively, we used a set of mathematical equations to process the raster images using a raster image calculator (Raster calculation tool within ArcMap 10.8). Landsat band 4, band 5, and thermal band 10 were used. Band 4 and band 5 used for calculating the Normalized Difference Vegetation Index (NDVI). There are six steps for estimating Landsat 8 LST [25,26,27] and they are described below.

(1)

The calculation of top-of-atmosphere radiance (TOA): The digital number (DN) of each pixel in the original Landsat 8 images was used to calculate the TOA radiance of the corresponding pixel according to the following formula provided by USGS:

$$TO{A}_{\lambda }=M_L\times Q_{cal}+A_L$$

(1)

where $TO{A}_{\lambda }$ represents the TOA (W m⁻² srad $\mathsf{\mu }\mathrm{m}$); $M_L$ is the band 10-specific multiplicative rescaling factor; $Q_{cal}$ is the DN values of band 10; and $A_L$ corresponds to the band-specific additive rescaling factor. The values for the above required parameters are available from the metadata file.

(2)

The conversion of TOA to brightness temperature (BT): The TOA radiance was converted to BT by Equation (2):

$$BT=K_2/\left(ln\left(\left(K_1/TOA\right)+1\right)\right)-273.15$$

(2)

where $K_2$ and $K_1$ represent the specific band thermal conversion constants, respectively. To obtain results in degrees Celsius, the radiation temperature needs to be adjusted by adding absolute zero (about—273.15 °C).

(3)

The NDVI calculation: The calculation of NDVI is important because the vegetation proportion ($P_{\nu }$) is highly related to NDVI and the emissivity (ε) that is related to $P_{\nu }$ can be calculated.

$$NDVI=\frac{\left(b_5-b_4\right)}{\left(b_5+b_4\right)},$$

(3)

where NDVI is the Normalized Difference Vegetation Index and b₅ and b₄ represent the fifth band (0.845–0.885 μm) and the fourth band (0.630–0.680 μm) of the Landsat 8 images, respectively.

(4)

The vegetation proportion calculation:

$$P_{\nu }=square\left(\frac{\left(NDVI-NDV{I}_{min}\right)}{\left(NDV{I}_{max}-NDV{I}_{min}\right)}\right)$$

(4)

where $NDV{I}_{min}$ is the minimum value of NDVI and $NDV{I}_{max}$ is the maximum value of NDVI.

(5)

Emissivity calculation:

$$\epsilon =0.004\times P_{\nu }+0.986$$

(5)

(6)

LST estimation:

$$LST=\frac{BT}{\left(1+\left(\frac{0.0015\times BT}{1.4388}\right)\times ln\left(\epsilon \right)\right)}$$

(6)

All calculations were completed using the ArcGIS raster calculation tools.

2.3.2. Field Measurements and Validation

The sampling sites were ensured not to be covered by vegetation before measurement, and the chamber was required to be close to the ground to ensure accurate reading of carbon flux. We inserted the chamber into the surface soil to a depth of 3–5 cm. The instrument recording time for each chamber was adjusted to 1 s, and the measurement time interval was 2–5 min, repeated five times. Due to the difference in daily variation of FSCO₂, the measurements between 10:00 AM to 3:00 PM (local time) were selected and used as the field monitoring FSCO₂ data in this study. The final FSCO₂ calculation equation [28] is as follows:

$$FSC{O}_2=\rho \times \frac{V}{S}\times \frac{\Delta C}{\Delta T}$$

(7)

where the FSCO₂ is the carbon release between the land surface and air in forest soil (g C·m⁻²·day⁻¹); $\rho$ is the density of CO₂; S and V represent the area (m²) and the volume (m³) of the chamber, respectively; ∆C is the changing trend of CO₂ concentration in soil with time; and ∆T is the time interval. We randomly selected 70% of FSCO₂ data for model operation, and the remaining 30% were used for model validation.

2.3.3. STI-FM Fusion Model and Validation

The STI-FM model is a fusion method proposed by Khaled hazaymeh in 2015 [17]. The aim of the model is to use two different sources of satellite remote sensing images to generate high-quality spatiotemporal resolution fusion products. The model is based on two assumptions. First, there is a linear relationship between two inconsistent MODIS LST images; and second, the LSTs obtained from Landsat 8 and MODIS images at a specific time (such as $T_1$ or $T_2$) are similar. The core idea of the model is to use the linear relationship between the two MODIS LST products to generate Landsat 8 LST prediction data of high time series. In other words, the model uses the linear relationship between the MODIS LST in $T_1$ and $T_2$ to generate a synthetic Landsat 8 LST image in $T_2$ using the Landsat 8 LST image in $T_1$.

$$MODIS\left(T_2\right)=a\times MODIS\left(T_1\right)+c$$

(8)

$$Synth\_L8\left(T_2\right)=a\times L8\left(T_1\right)+c$$

(9)

where $MODIS\left(T_2\right)$ and $MODIS\left(T_1\right)$ are two consecutive MODIS LST images; a and c are the slope and intercept between the $MODIS\left(T_2\right)$ and $MODIS\left(T_1\right)$, respectively; $L8\left(T_1\right)$ is the Landsat 8 image collected from the same time with $MODIS\left(T_1\right)$; and $Synth\_L8\left(T_2\right)$ is the fusion image in $T_2$ at the same time and spatial resolution as $L8\left(T_1\right)$.

Hazaymeh applied the STI-FM model to the semi-arid area of the Middle East and Jordan [17]. When applying it to the subtropical monsoon climate region, it is necessary to verify the applicability of the model. Previously, scholars applied the model to the Great Bay area of Guangdong, Hong Kong, and Macao and verified its feasibility [16]. The model was performed in the Google Earth Engine (GEE) platform, using R language and ArcGIS 10.8 software. The accuracy of the model is evaluated by assessing the adjusted $R^2$ and root mean square error ($RMSE$) using R language.

$$R^2={\left[\frac{{\displaystyle \sum }\left(L_{\left(t\right)}-{\overline{L}}_{\left(t\right)}\right)\left(S_{\left(t\right)}-{\overline{S}}_{\left(t\right)}\right)}{\sqrt{{\left({\displaystyle \sum }{L}_{\left(t\right)}-{\overline{L}}_{\left(t\right)}\right)}^2}\sqrt{{\left({\displaystyle \sum }{S}_{\left(t\right)}-{\overline{S}}_{\left(t\right)}\right)}^2}}\right]}^2$$

(10)

$$RMSE=\sqrt{\frac{{\displaystyle \sum }{\left[S_{\left(t\right)}-L_{\left(t\right)}\right]}^2}{n}}$$

(11)

where $L_{\left(t\right)}$ is the actual Landsat 8 LST value; $S_{\left(t\right)}$ is the Landsat 8 LST value synthesized by the model; ${\overline{L}}_{\left(t\right)}$ is the average value of the actual Landsat 8 LST; ${\overline{S}}_{\left(t\right)}$ is the average value of Landsat 8 LST synthesized by the model; and N = 10,000.

2.3.4. FSCO₂ Estimation and Validation

The estimation of FSCO₂ based on remote sensing data involves the following three steps: Firstly, a normality test is conducted on the field-measured FSCO₂ data and the fused LST product test to explore seasonal, monthly, and daily changes of FSCO₂. Then, the constrained maximum likelihood method is applied to establish three linear mixed effect models under different time scales (season, month, and day). The linear mixed models include the null model, where the random intercept as a fixed predictor is grouped according to the level of daily and monthly factors, and the Daily random effect model, where the fusion of remote-sensing-derived LST is used as a fixed predictor and the daily factor level is used as a random intercept effect. The fusion remote sensing LST as a fixed predictor plus the daily factor level provides a random intercept effect; the composition of the Monthly random effect model and Seasonally random effect models are similar to that of the Daily random effect model, but the monthly and seasonal factor levels are used as the random intercept effects, respectively, in the two models. Lastly, the three models were compared with their corresponding null models utilizing the Akaike information criterion (AIC) and Bayesian information criterion (BIC) to select the best optimum model.

Marginal R² (mR²) and conditional R² (cR²) were used to evaluate these models in estimating seasonal, monthly, and daily FSCO₂ [29]. Marginal R² quantifies the proportion of variance explained by the fixed prediction of Landsat 8 LST, while condition R² evaluates the proportion of variance captured by the fixed prediction of Landsat 8 LST and the seasonal, monthly, and daily factor levels (random factors). The accuracy and efficiency of the Daily random effect model, Monthly random effect model, and Seasonally random effect model were evaluated using the validation data set. The accuracy of these models was evaluated by assessing their adjusted R² and RMSE.

3. Results

3.1. Field-Measured FSCO₂ and Related Environmental Variables Analysis

Figure 3 shows the seasonal and daily changes in FSCO₂ in the CHD and the GY. It is evident that there are significant seasonal differences in FSCO₂ at both study sites, with higher values in autumn than in winter (Figure 3 CHD Season and GY Season). Additionally, the average FSCO₂ values in Site site GY are smaller than that in site CHD, indicating that FSCO₂ in the study area varies both temporally and spatially. When abnormal values are excluded, the FSCO₂ in GY during the dry season ranges from 0.15 to 39.36 g C m⁻² day⁻¹, with an average of 7.83 g C m⁻² day⁻¹. In contrast, the FSCO₂ in CHD during the dry season ranges from 0.97 to 41.09 g C m⁻² day⁻¹, with an average of 12.73 g C m⁻² day⁻¹. This indicates that the FSCO₂ values in the study sites are spatially heterogeneous, depending on time and space.

Figure 4 shows the relationship between the field-measured FSCO₂ and the soil respiration, topsoil temperature, air temperature, and relative humidity. Results show that the correlation coefficients between the FSCO₂ and surface air temperature in the GY and CHD sites are 0.46 (p < 0.001) and 0.70 (p < 0.001), respectively. The topsoil temperature in the GY and CHD sites has stronger positive correlation with the FSCO₂ (R² = 0.72 and R² = 0.70) at the 0.001 significance level. FSCO₂ is also significantly correlated with the relative humidity in GY; the correlation coefficient is 0.4. The results also show that changes in soil respiration do not have a noticeable impact on the level of FSCO₂, and that soil temperature and near-surface air temperature are crucial environmental factors affecting the FSCO₂ in the subtropical forest during the dry season.

3.2. Multisource Remote Sensing LST Fusion

3.2.1. Fusion LST Datasets and Accuracy Assessment

Figure 5A shows a qualitative comparison of Landsat 8 LST and Fusion LST images on 28 September 2019, highlighting their similarity in features. The study also investigated LST at three different elevations during the study period and found that the Fusion LST image accurately predicted LST in variable topographic conditions. Histograms were generated for actual Landsat 8 LST and Fusion LST images for the entire study area and revealed their similarities (Figure 5B). For quantitative evaluation, Fusion LST and MYD LST of the entire study area were plotted on 28 September 2019, and further Fusion LST and topsoil temperature in sampling chamber points were plotted during the study period (Figure 5C,D). Strong relationships were found between the variables of interest, with R² values of 0.77 and 0.60, and RMSE values of 0.33 and 2.24, respectively. Additionally, the close relationship between the regression line of MYD-LST and Fusion LST and the 1:1 line indicates a strong correlation between the two datasets.

3.2.2. The Spatiotemporal Variations of LST in Landsat 8 LST, MODIS LST and Fusion LST

Spatial variations of LST show that the distribution of FSCO₂ has significant spatial patterns in all three multisource remote sensing LST images in different study periods and that surface topography is the major factor responsible for the spatial variations of FSCO₂: as the altitude rises, the LST decreases. The LST variations also significantly differ during different seasons of the year. The monthly mean temperatures of September, October, and December during the dry season are 25 °C, 20 °C, and 15 °C in the study region (Figure S1), respectively. Moreover, compared with the MODIS LST (500 m), the Fusion LST (100 m) has significantly higher spatial resolution, which improved markedly the FSCO₂ estimation at large regional scale. The accuracy of fusion LST data satisfied the conditions for simulating FSCO₂.

Figure 6 shows the spatiotemporal changing patterns of LST in autumn, winter, and the dry seasons as detected by Landsat and MODIS images and Fusion LST products of the study area. The figure illustrates that there is considerable variation in the accuracy of slope when using multisource satellite remote sensing products. The changing trends of LST at the pixel level obtained using Landsat images have higher spatial resolution, but the trend analysis is not accurate due to the low time resolution. This suggests that using a combination of satellite imagery from different sources can improve the accuracy of LST mapping, but the low time resolution of the images may still affect the accuracy of the results.

Furthermore, Figure 6 shows the trends of LST in the subtropical forests in different seasons. The changing slopes of LST in winter are larger than in autumn, which is likely due to the significant decrease in LST in subtropical forests during the winter season. The spatial distribution of the trend variance in LST reveals distinct patterns over different periods. Specifically, the results indicate that LST in high-altitude areas differs more significantly across space compared with low-altitude areas. This is likely because the forest canopy in high-altitude areas is denser than in low-altitude areas, resulting in LST being more similar to near-surface air temperature. As a result, the changing dynamics of LST in high-altitude areas are more significant than in low-altitudinal areas.

3.3. FSCO₂ Simulation

3.3.1. The Construction of the Linear Mixed Based FSCO₂ Inversion Model

The similarities in variance between the Fusion LST variable and residual variance across different months indicate that the changes in Fusion LST are mainly caused by random factors rather than by other fixed factors (Table S1). The degree of change in Fusion LST across different months is relatively small (0.04) in the observed data, and this difference is mainly due to random error, suggesting that the variations in Fusion LST between months are mainly caused by random factors such as measurement error. Additionally, the differences in Fusion LST between months account for 46.8% of the total variability in the data, and the remaining variability can be explained by random errors in the model.

For the fixed effect, the slope parameter of 0.21 and the intercept of 0.54 suggest that there is a linear relationship between the independent and dependent variables. The data indicate that there is a significant difference (p < 0.001) in mean Fusion FSCO₂ between measurements taken from September to January in the dry season. Additionally, the variance in Fusion LST within months is relatively low (0.03) compared with the variance between autumn and winter (0.05). This suggests that there is less variation in Fusion LST within each month than between autumn and winter. It is noticeable that the variation of Fusion LST by month accounted for 34.5% of the error left in the model, indicating that when the predictor variable (Fusion LST) is grouped by autumn and winter, there is no significant difference in FSCO₂ measurement between autumn and winter. This also suggests that the Fusion LST model does not vary significantly between these two seasons. It means that the model is effective in explaining the variance in Fusion FSCO₂ between months, but less effective in explaining the variance between autumn and winter (Table S2).

The variance in Fusion LST between days was high (variance = 0.07) compared with that within months (variance = 0.007). This suggests that there is more variation in Fusion LST between days than within each month. The variation of Fusion LST by day accounted for 9.72% of the error left in the model. For the fixed effect, the slope parameter was 0.13 while the intercept was 0.48. This suggests that there is a linear relationship between the independent and dependent variables, but the slope of 0.13 is relatively small, which means that the change in dependent variable is not large for one unit change in the independent variable. The intercept is 0.48, which means when the independent variable is zero, the dependent variable has a value of 0.48. The result also states that there is no significant difference (p < 0.42) in mean Fusion FSCO₂ measurements. This means there is not enough evidence to show that a meaningful difference in the mean Fusion FSCO₂ measurements between different groups or conditions exists. In other words, it suggests that the model is not effective in explaining the variance in Fusion FSCO₂ measurements between days and within months, and that the data do not provide enough evidence for a meaningful difference in the mean Fusion FSCO₂ measurements between different groups or conditions (Table S3).

3.3.2. Model Validations

Table 2. Evaluation of linear mixed effect models, outlined by random intercept and Fusion LST segmented by day (between 28 September 2019 to 4 January 2020), month (September, October, November, December, and January), and season (Autumn, Winter and Dry) of the year. Models were assessed using the Akaike information criterion (AIC), Bayesian information criterion (BIC), the p-value of a chi-squared test, and both marginal and conditional R-squared (mR² and cR²).

Model Type	Model Description	AIC	BIC	p-Value	mR2	cR2
None random effect	FSCO2~1 + (1|day)	-	-	-
	FSCO2~1+ (1|month) FSCO2~1 + (1|season)	-	-	-
Daily random effect	FSCO2~LST + (1|day)	220.75	208.99		0.19	0.904
Monthly random effect	FSCO2~LST + (1|month)	−34.76	−23.00	<0.001	0.47	0.553
Seasonally random effect	FSCO2~LST + (1|season)	984.47	996.24	<0.217	0.07	0.353

shows the performance of the Daily random effect model, Monthly random effect model, and Seasonally random effect model in comparison with their respective null models and in comparison with each other. The results indicate that the Daily random effect model and Monthly random effect model comparatively have lower AIC and BIC values than the Seasonally random effect model. This suggests that the Daily and Monthly models have better explanatory power and are more parsimonious than the seasonal model. Meanwhile, the Daily random effect model and Monthly random effect model were significantly improved over the intercept-only model.

The results of a statistical analysis of three different models indicate that the use of fixed predictors (such as Daily, Monthly, and Seasonally) can explain variations in the dependent variable (FSCO₂) (Table 2). For example, the three models performed differently in terms of how well their fixed predictors (mR²) or combination of fixed predictors and random factors (cR²) explained the variations in FSCO₂. The fixed predictor of the Daily random effect model captured 19% of the variability in FSCO₂, the Monthly random effect model explained 47%, and the Seasonally random effect model explained 7%. However, when the proportions of variance captured by fixed predictors and random factors were combined, the Daily random effect model performed better than the Monthly and Seasonally random effect models. The former could explain 90.4% of the total variance in the model, while the Monthly random effect model and Seasonally random effect model explained only 55.3% and 35.3%, respectively. It is noteworthy that the Daily random effect model and Seasonally random effect model do not deviate from the assumption of common variance of linear regression models. A likelihood-ratio test showed that there was no statistically significant difference between a heteroscedastic model and each of the Monthly random effect model, Daily random effect model, and Seasonally random effect model.

A detailed analysis of the results showed that the discrepancies between the estimated values and the actual observed FSCO₂ varied depending on which model was used (Figure 7 and Figure 8). This indicates that the three models had different levels of accuracy in explaining the variations of the dependent variable of FSCO₂. Before FSCO₂ reached 0.5 g C m⁻²/day⁻¹, the Daily random effect model overestimated FSCO₂. When FSCO₂ reached 0.75 g C m⁻²/day⁻¹ in the Monthly random effect model, FSCO₂ values were estimated properly in the Seasonally random effect model. Consequently, the Daily random effect model produces an RMSE of 0.22, while the RMSE is 0.28 in the Monthly random effect model and 0.9 in the Seasonally random effect model. Although the Daily random effect model and the Monthly random effect model exhibited similar RMSE, the Daily random effect model demonstrated a stronger correlation between the estimated and observed FSCO₂ than the Monthly random effect model (adjusted R² = 0.44, 0.9, and 0.19 for the Daily random effect model, Monthly random effect model, and Seasonally random effect model, respectively) (Figure 7). The Monthly random effect model that performed optimally had an explanation rate of 55.3% (conditional R² = 0.55 and t value > 1.9) for FSCO₂ variability and yielded the smallest deviation from observed FSCO₂. While all three models have biases, the Taylor diagram for the monthly model is closer to the observed data points (Figure 8D). The quantiles of both the monthly and daily models have lower errors (Figure 8A–C), which indicates that they have strong predictive performance. Overall, the results suggest that the Monthly random effect model is the best model and can explain the most variation in the dependent variable of FSCO₂. It also adheres to the assumption that the variance of the errors is constant across all levels of the independent variable in a linear regression model.

4. Discussion

4.1. Spatial Heterogeneity of FSCO₂ Variations and Relationship between the Abiotic Factors

Based on the linear mixed effect models (e.g., Monthly random effect model), we mapped FSCO₂ spatial distribution during the study periods, and further calculated the pixel-based mean and slope values of FSCO₂ in the dry, autumn, and winter seasons to understand the spatiotemporal distribution characteristics of FSCO₂ in the subtropical region during the dry season (Figure 9). The mean FSCO₂ values in different seasons have small FSCO₂ in the north and center of the study areas compared with other areas. However, the mean FSCO₂ is variable in different periods. The highest FSCO₂ occurred in Autumn and Winter. The variation trends of FSCO₂ also have significant variation in different study periods, and the changing dynamics is the highest in autumn compared to other periods.

The predicted FSCO₂ and abiotic factors have significant relationships in dry seasons. Moreover, the degree of relevance is different in GY and CHD. Figure 10 shows that predicted FSCO₂ and topsoil temperature and near-surface air temperature have significant positive relationships; however, the relationship between the FSCO₂ and soil respiration is very weak in the GY site, while no significant relationship exists in the CHD site. Moreover, the correlation coefficients (p < 0.001) in the GY site were significantly higher than those in the CHD, illustrating that topography influences the relationship between the FSCO₂ and abiotic environmental factors.

4.2. Spatiotemporal Dynamics of Soil CO₂ Efflux of Subtropical Forest in Dry Seasons

According to the results, the FSCO₂ during the dry season may offset a significant amount of CO₂ assimilated during the growing season, potentially accounting for 3–50% of annual carbon emissions. Therefore, FSCO₂ in the dry season is crucial to determine the annual carbon cycle [4,22,30,31]. In our study, the subtropical forest in dry seasons was targeted for using the remote sensing fusion method to acquire high-quality spatiotemporal products of LST and to estimate FSCO₂ in dry seasons. Our results indicated that FSCO₂ trended downward from 28 September to 4 January of the following year, which is consistent with the results of Chen et al. [22], who concluded that FSCO₂ has distinct variation during the dry season. Our results further revealed the spatiotemporal variability of FSCO₂ in dry seasons at high spatiotemporal resolutions, and significantly improved the distribution accuracy of FSCO₂ estimation. For example, at the temporal scale, we estimated daily FSCO₂ in the LXH reservoir basin. At spatial scales, we estimated 100 m spatial resolution of FSCO₂ while significantly decreasing mixed-pixel uncertainty that would arise from coarse-spatial-resolution images. In addition, our estimated mean maximum and minimum values of FSCO₂ in the autumn, winter, and dry seasons are lower than the estimated values by Chen et al. [22]. This discrepancy could be due to the following two reasons. First, there is a difference in the data sources for FSCO₂ estimation. Chen et al. [22] used MODIS images at coarse 500 m spatial resolution that may have contained many mixed pixels that biased the true estimation of FSCO₂. Second, Chen et al. [22] assumed that the FSCO₂ estimation values were affected by plant productivity and soil moisture but not land surface temperature. In our study, we aimed to mainly use the land surface temperature to improve the inversion accuracy of FSCO₂ estimation at high temporal and spatial resolutions at regional scales. Nevertheless, our future study will add other biotic and abiotic factors, such as soil moisture and plant productivity, to further improve the model’s accuracy. Moreover, our results revealed that the long spatiotemporal dynamics of FSCO₂ experience daily, monthly and seasonal variability.

4.3. Remote Sensing Based Soil CO₂ Efflux Inversion

Remote sensing products are of great potential for the estimation of FSCO₂ emission at long-time and large spatial scales [32]. Earth observation technology allows for the collection of data over large areas and over long time periods, making it well-suited for monitoring and estimating emissions at regional or global scales. Additionally, the use of remote sensing data can help to reduce the cost and logistic challenges associated with ground-based measurements. Freely accessible remote sensing images like Landsat, MODIS, and ASTER were commonly used in previous studies for estimating FSCO₂ [10,32,33] and obtained a relatively optimistic accuracy in LST (

Table 3. Previous studies of remote-sensing-based soil CO₂ efflux monitoring during 2000–2022.

Authors	Sites	Method	Satellite Data and Resolutions	Inversion Spatial Scales	Highest Inversion Accuracy	Published Year
Kimball et al. [34]	Tundra forest	terrestrial carbon flux (TCF) model	MODIS and AMSR-E	Continent scale	0.89	2009
Huang et al. [35]	Broadleaf forest site (Midwest USA)	Statistical model	MODIS LST/500 m	Basin scale		2014
Wu et al. [3]	Canadian boreal black spruce stand	Linear regression model	MODIS LST & NDVI/500 m	Landscape scale	0.78	2014
Huang et al. [11]	FLUXNET forest	Remote-sensing-based model	MODIS LST/500 m	Site scale		2015
Huang et al. [36]	Croplands	support vector regression	Landsat 8 images	County level	0.73	2017
Ben Bond-Lamberty [37]		Artificial neural network model		Global scale		2018
Crabbe et al. [32]	Forest	Linear mixed model	Landsat 8 LST/30 m	Patch scale		2019
Warner et al. [38]	Bamboo forest	Quantile-based digital soil mapping	DEM/2 m	Basin scale	0.64	2019
Huang et al. [39]	Global	biome-specific statistical model	MODIS/1 km	Global scale		2020
Xu et al. [40]	Forest	Improved downscale model	MODIS, Landsat 8 OLI/TIRS	Regional scale	0.47	2020
Chen et al. [22]	Tropical forest	Random forest	MODIS/500 m	Basin scale	0.88	2021
Burdun et al. [33]	Peatlands	Model	Landsat/MODIS LST/1 km	Regional scale	0.67	2021

). Wu et al. [3] reported long-time-series measurements of FSCO₂ using the relationship between the FSCO₂ and MODIS LST. Richard et al. [32] used nighttime and daytime MODIS LST data to monitor the FSCO₂ for shorter time periods. Satellites equipped with thermal infrared detectors can be a valuable tool for measuring soil temperature at a large spatial scale, particularly in forested areas. This method can provide more accurate and extensive understanding of land surface temperature, and its potential impact on the environment, by circumventing the limitations of traditional in situ methods. It provides a more efficient and comprehensive approach to monitor soil temperature at landscape or ecosystem scales [3]. Richard et al. [32] assumed that satellite-remote-sensing-based LST values obtained from monitoring of the forest canopy can reflect forest canopy temperature. The forest canopy size plays a crucial role in determining the accuracy of LST estimation, as the heat retention capacity of the dense canopy during peak growth seasons is contingent on the air conditions in its immediate vicinity. Therefore, with an increase in forest canopy area, LST is close to the air temperature [32]. In our study, the forest canopy of the entire study area was greater than 98%. The Landsat-MODIS LST inversion therefore has a relatively positive relationship with air temperature, as shown in Figure 10.

Multisource remote sensing image fusion technology improves the quality of data by effectively fusing data from different sources and features, and by leveraging the unique advantages of various remote sensing data in terms of both spatial and temporal resolution. The calibration and verification of the fusion results can usually be achieved by using ground observation data. When it comes to using satellite technology to measure temperature on the earth’s surface, we can obtain more accurate results by combining information from several different sources. We can achieve better precision by analyzing images that were acquired at different points in time, from varying perspectives, and at varying levels of radiation intensity [32]. In this study, that is why we used both Landsat and MODIS LST products to estimate high spatiotemporal FSCO₂ at basin scales.

4.4. Accuracy Analysis of Remote Sensing Modeling for FSCO₂ Estimation

The high spatiotemporal simulation of the carbon flux between the surface and the atmosphere depends heavily on the field in situ observation of carbon flux. However, regional or general climate models ignore scale issues when monitoring carbon fluxes at a coarser scale. This limits our large-scale and high spatiotemporal retrieval of carbon flux at the soil-air interface. Our method provides an opportunity to estimate large-scale FSCO₂ from surface temperatures determined from the fusion of multisource satellite imagery.

The accuracy of FSCO₂ estimation based on satellite remote sensing images is subject to several factors, including satellite access time, specific distribution characteristics, regional topographic conditions, vegetation coverage, soil and air temperature, etc. [32]. The assessment of soil FSCO₂ in this study shows that soil moisture and temperature in the subtropical forests in the dry season are the main drivers of FSCO₂ change. It also indicated that the relationship between soil moisture and FSCO₂ is more complex than the relationship between soil temperature and FSCO₂. In other words, different degrees of soil wetness may apply different effects on the FSCO₂ efflux. The soil moisture (<50%) has a positive effect on FSCO₂ efflux for the subtropical forests in the dry season [22]. Another study found that the accuracy of estimating vegetation indices such as NDVI, LAI, and LST from satellite imagery is subject to soil moisture. The models developed in this study took this into account [3]. Previous studies have shown that oversaturated soil has a limiting effect on soil CO₂ release [3,32,33]. It is generally recognized that soil moisture in subtropical humid regions is relatively high and rainfall is sufficient. Therefore, for this reason, we have chosen the dry season (September to January of the following year) as the observation period. There was no rain in the study area during the study period, and the average soil moisture of the two sites was 18.10% and 20.99%, both less than 50%. Therefore, there was a positive correlation between FSCO₂ and soil moisture. Specifically incorporating soil moisture as a model parameter into a monthly model (Monthly random effect model) (at least when considering summer FSCO₂ variability) is necessary. However, the lack of high-quality soil moisture data with high spatial resolution limited this study. Therefore, improving the model using soil moisture from satellite remote sensing data is needed in future research.

4.5. Limitations and Future Works

This work attempts to deal with the accurate calculation of FSCO₂ at large spatial scales in a subtropical forest based on high spatiotemporal LST inversions. In other words, it provides an opportunity to estimate FSCO₂ at large spatial scales based on the surface temperature as determined by satellite remote sensing images. Although remote sensing fusion technology enables high temporal resolution for FSCO₂ estimation, there are still several limitations that need further addressing. Firstly, field measurement of FSCO₂ used chambers at five different times a day. LST is very sensitive to the air and soil temperature, and although we acquired quite satisfactory results from models, a better method should use the data with timing corresponding well with the satellite accessing time. Secondly, this paper collected observation data from two sites over a period of five months. The limited duration of the observations may have resulted in a relatively low accuracy of the model. Future research will attempt to improve the model by expanding the spatial and temporal scale of the observation data.

5. Conclusions

Utilizing satellite remote sensing images for FSCO₂ estimation has significance for accurate calculation of the carbon budget balance. Accurate spatiotemporal FSCO₂ estimation at regional scales is an urgent issue that needs immediate solutions for regional climate models. In this research, we estimated soil CO₂ efflux using a combination of Landsat and MODIS imagery and constructed linear mixed effect models that account for daily, monthly, and seasonal variability in Landsat-MODIS Fusion LST in the dry season. By using daily, monthly, and seasonal FSCO₂ measurements as separate random factors in the linear mixed effect model, three different FSCO₂ models were built. Our findings revealed that the monthly random effect model, which utilizes Landsat 8 LST and MODIS LST fusion images, can best explain the forest FSCO₂ dynamics at a regional scale. There is a strong positive correlation between the predicted FSCO₂ and abiotic environmental factors such as air and soil temperature. Future research can enhance these models by incorporating other biotic and abiotic factors such as plant productivity, soil moisture, etc. The estimation of soil CO₂ emissions from subtropical forests based on remote sensing is still in its infancy, and, to a large extent, this study deepened the knowledge of this field.