Abstract: Objective　Satellite remote sensing has emerged as a pivotal tool for capturing large-scale, high-resolution atmospheric carbon dioxide (CO2) concentration data, which is indispensable for tracking the progress of China's carbon neutrality goals. As global climate change intensifies, accurate and continuous monitoring of regional CO2 dynamics has become increasingly critical for formulating targeted mitigation strategies. However, the acquisition of comprehensive, high-quality data over large regions within short timeframes (e. g., one month or one year) is still constrained by various factors, including the limited number of CO2-monitoring satellites, narrow swath widths, persistent cloud cover, unfavorable lighting conditions, and aerosol scattering effects. The Beijing-Tianjin-Hebei (BTH) region, one of the most economically advanced and densely populated areas in China, is characterized by intensive industrial activities, high energy consumption, and carbon emissions accounting for approximately 11% of the national total. In addition, its unique geographical location, complex terrain, and diverse land-use patterns further amplify the complexity of CO2 spatiotemporal distributions. Therefore, a comprehensive understanding of the spatiotemporal dynamics of CO2 column concentration (XCO2) in this region is critical for formulating effective strategies to achieve national carbon peak and carbon neutrality targets. Methods　To address data gaps and generate spatially continuous XCO2 data, this study employed the empirical Bayesian 265 Kriging interpolation method to fill missing values in OCO-2 satellite XCO2 datasets. Monthly XCO2 data with a high spatial resolution of 0.1° (approximately 11 km) for the BTH region from 2015 to 2021 were produced, enabling fine-scale analysis of regional variations. The primary data source was the OCO-2 Level 2 XCO2 product (OCO2_L2_Lite_FP v10r), renowned for its high accuracy and spatial detail, supplemented by multi-source auxiliary data including the normalized difference vegetation index (NDVI) from the MOD13C2 product, precipitation and temperature data from the ERA5 reanalysis dataset, and anthropogenic carbon emissions data from the Open-source Data Inventory for Anthropogenic CO2 (ODIAC). Validation of the interpolated data was rigorously conducted using ground-based XCO2 measurements from Xianghe station of the Total Carbon Column Observing Network (TCCON), a benchmark for high-precision atmospheric CO2 monitoring. Additionally, time series decomposition and correlation analysis were performed to systematically explore the spatiotemporal characteristics of XCO2 and its complex relationships with environmental factors (vegetation, precipitation, temperature) and anthropogenic drivers (fossil fuel emissions). Results and Discussion　 The results demonstrated that the empirical Bayesian Kriging interpolation significantly enhanced data coverage, with average increases of 95.4%, 88.36%, and 69.19% on monthly, seasonal, and annual scales, respectively, effectively addressing the data sparsity issue inherent in satellite observations. The validation results indicated strong consistency between the interpolated data and the ground observations, with a determination coefficient (R²) of 0.84, a root mean square error (RMSE) of 1.91 μmol/mol, and a mean absolute error (MAE) of 1.40 μmol/mol, all of which were superior to those of the original OCO-2 dataset (R² = 0.79, RMSE = 2.20 μmol/mol, MAE = 2.05 μmol/mol). Temporally, the XCO2 concentrations in the BTH region exhibited a continuous annual growth trend, increasing from 400.34 μmol/mol in 2015 to 416.67 μmol/mol in 2021, with an average annual increase of 2.33 μmol/mol, reflecting the cumulative impact of anthropogenic emissions despite regional emission reduction efforts. Seasonally, the XCO2 concentrations in summer and autumn were distinctly lower than those in spring and winter, peaking in April and reaching their minimum in August, which was consistent with the seasonal rhythm of vegetation photosynthesis and energy consumption. Spatially, the Taihang- Yanshan Mountains served as a natural boundary, resulting in that the XCO2 values were lower in the northern mountainous areas (dominated by forests and grasslands with strong carbon sink capacity) and higher in the southern plains (characterized by intensive industrial and agricultural activities). Correlation analysis revealed that NDVI (−0.37), precipitation (−0.33), and temperature (−0.36) were negatively correlated with monthly average XCO2, with even stronger correlation with the seasonal component of XCO2 (NDVI: −0.77; temperature: −0.57; precipitation: −0.62), highlighting the critical role of vegetation in regulating seasonal CO2 fluctuations. Anthropogenic carbon emissions showed a positive correlation with the monthly XCO2 (0.45) but a negative correlation with the seasonal component ( − 0.35), primarily due to that the enhanced vegetation photosynthesis efficiently absorbed CO2 during the vegetation growing season (April-August) despite the increased emissions from human activities. Conclusions　This study confirms the effectiveness of the empirical Bayesian Kriging method in generating high-precision, spatially continuous XCO2 data for complex regions like BTH. The identified spatiotemporal patterns of XCO2 and their driving factors provide valuable scientific support for regional carbon emission reduction policies and the achievement of carbon peak and carbon neutrality targets. The findings underscore the need for formulating differentiated low-carbon strategies based on regional characteristics, such as strengthening vegetation carbon sinks in mountainous areas through ecological restoration and optimizing industrial emissions in plain regions via technological innovation and energy structure adjustment. These insights can provide guidance for targeted and efficient climate action, contributing to the sustainable development of the BTH region and the broader national carbon neutrality agenda.

Key words: OCO-2, CO2 column concentration, spatiotemporal distribution, Kriging interpolation